%DECK DSS004
      function [ux]=dss004(xl,xu,n,u)
%...
%...  FUNCTION DSS004 COMPUTES THE FIRST DERIVATIVE, U , OF A
%...                                                  X
%...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL
%...  FIVE-POINT, FOURTH-ORDER FINITE DIFFERENCE APPROXIMATIONS
%...
%...  ARGUMENT LIST
%...
%...     XL      LOWER BOUNDARY VALUE OF X (INPUT)
%...
%...     XU      UPPER BOUNDARY VALUE OF X (INPUT)
%...
%...     N       NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE
%...             BOUNDARY POINTS (INPUT)
%...
%...     U       ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT
%...             THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS
%...             TO BE COMPUTED (INPUT)
%...
%...     UX      ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL
%...             VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS
%...             (OUTPUT)
%...
%...  THE MATHEMATICAL DETAILS OF THE FOLLOWING TAYLOR SERIES (OR
%...  POLYNOMIALS) ARE GIVEN IN FUNCTION DSS002.
%...
%...  FIVE-POINT FORMULAS
%...
%...     (1)  LEFT END, POINT I = 1
%...
%...                                   2            3            4
%...  A(U2 = U1 + U1  ( DX) + U1  ( DX)  + U1  ( DX)  + U1  ( DX)
%...                X   1F      2X  2F       3X  3F       4X  4F
%...
%...                      5            6            7
%...           + U1  ( DX)  + U1  ( DX)  + U1  ( DX)  + ...)
%...               5X  5F       6X  6F       7X  7F
%...
%...                                   2            3            4
%...  B(U3 = U1 + U1  (2DX) + U1  (2DX)  + U1  (2DX)  + U1  (2DX)
%...                X   1F      2X  2F       3X  3F       4X  4F
%...
%...                      5            6            7
%...           + U1  (2DX)  + U1  (2DX)  + U1  (2DX)  + ...)
%...               5X  5F       6X  6F       7X  7F
%...
%...                                   2            3            4
%...  C(U4 = U1 + U1  (3DX) + U1  (3DX)  + U1  (3DX)  + U1  (3DX)
%...                X   1F      2X  2F       3X  3F       4X  4F
%...
%...                      5            6            7
%...           + U1  (3DX)  + U1  (3DX)  + U1  (3DX)  + ...)
%...               5X  5F       6X  6F       7X  7F
%...
%...                                   2            3            4
%...  D(U5 = U1 + U1  (4DX) + U1  (4DX)  + U1  (4DX)  + U1  (4DX)
%...                X   1F      2X  2F       3X  3F       4X  4F
%...
%...                      5            6            7
%...           + U1  (4DX)  + U1  (4DX)  + U1  (4DX)  + ...)
%...               5X  5F       6X  6F       7X  7F
%...
%...  CONSTANTS A, B, C AND D ARE SELECTED SO THAT THE COEFFICIENTS
%...  OF THE U1  TERMS SUM TO ONE AND THE COEFFICIENTS OF THE U1  ,
%...           X                                                2X
%...  U1   AND U1   TERMS SUM TO ZERO
%...    3X       4X
%...
%...  A +   2B +   3C +   4D = 1
%...
%...  A +   4B +   9C +  16D = 0
%...
%...  A +   8B +  27C +  64D = 0
%...
%...  A +  16B +  81C + 256D = 0
%...
%...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
%...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
%...                                                            4X
%...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
%...               X
%...                                                         4
%...  U1  = (1/12DX)(-25U1 + 48U2 - 36U3 + 16U4 - 3U5) + O(DX )   (1)
%...    X
%...
%...     (2)  INTERIOR POINT, I = 2
%...
%...                                   2            3            4
%...  A(U1 = U2 + U2  (-DX) + U2  (-DX)  + U2  (-DX)  + U2  (-DX)
%...                X   1F      2X  2F       3X  3F       4X  4F
%...
%...                      5            6            7
%...           + U2  (-DX)  + U2  (-DX)  + U2  (-DX)  + ...)
%...               5X  5F       6X  6F       7X  7F
%...
%...                                   2            3            4
%...  B(U3 = U2 + U2  ( DX) + U2  ( DX)  + U2  ( DX)  + U2  ( DX)
%...                X   1F      2X  2F       3X  3F       4X  4F
%...
%...                      5            6            7
%...           + U2  ( DX)  + U2  ( DX)  + U2  ( DX)  + ...)
%...               5X  5F       6X  6F       7X  7F
%...
%...                                   2            3            4
%...  C(U4 = U2 + U2  (2DX) + U2  (2DX)  + U2  (2DX)  + U2  (2DX)
%...                X   1F      2X  2F       3X  3F       4X  4F
%...
%...                      5            6            7
%...           + U2  (2DX)  + U2  (2DX)  + U2  (2DX)  + ...)
%...               5X  5F       6X  6F       7X  7F
%...
%...                                   2            3            4
%...  D(U5 = U2 + U2  (3DX) + U2  (3DX)  + U2  (3DX)  + U2  (3DX)
%...                X   1F      2X  2F       3X  3F       4X  4F
%...
%...                      5            6            7
%...           + U2  (3DX)  + U2  (3DX)  + U2  (3DX)  + ...)
%...               5X  5F       6X  6F       7X  7F
%...
%...  -A +   B +  2C +  3D = 1
%...
%...   A +   B +  4C +  9D = 0
%...
%...  -A +   B +  8C + 27D = 0
%...
%...   A +   B + 16C + 81D = 0
%...
%...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
%...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
%...                                                            4X
%...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
%...               X
%...                                                        4
%...  U2  = (1/12DX)(-3U1 - 10U2 + 18U3 -  6U4 +  U5) + O(DX )    (2)
%...    X
%...
%...     (3)  INTERIOR POINT I, I NE 2, N-1
%...
%...                                        2             3
%...  A(UI-2 = UI + UI  (-2DX)  + UI  (-2DX)  + UI  (-2DX)
%...                  X    1F       2X   2F       3X   3F
%...
%...                          4             5             6
%...              + UI  (-2DX)  + UI  (-2DX)  + UI  (-2DX)  + ...)
%...                  4X   4F       5X   5F       6X   6F
%...
%...                                        2             3
%...  B(UI-1 = UI + UI  ( -DX)  + UI  ( -DX)  + UI  ( -DX)
%...                  X    1F       2X   2F       3X   3F
%...
%...                          4             5             6
%...              + UI  ( -DX)  + UI  ( -DX)  + UI  ( -DX)  + ...)
%...                  4X   4F       5X   5F       6X   6F
%...
%...                                        2             3
%...  C(UI+1 = UI + UI  (  DX)  + UI  (  DX)  + UI  (  DX)
%...                  X    1F       2X   2F       3X   3F
%...
%...                          4             5             6
%...              + UI  (  DX)  + UI  (  DX)  + UI  (  DX)  + ...)
%...                  4X   4F       5X   5F       6X   6F
%...
%...                                        2             3
%...  D(UI+2 = UI + UI  ( 2DX)  + UI  ( 2DX)  + UI  ( 2DX)
%...                  X    1F       2X   2F       3X   3F
%...
%...                          4             5             6
%...              + UI  ( 2DX)  + UI  ( 2DX)  + UI  ( 2DX)  + ...)
%...                  4X   4F       5X   5F       6X   6F
%...
%...   -2A -   B +   C +  2D = 1
%...
%...    4A +   B +   C +  4D = 0
%...
%...   -8A -   B +   C +  8D = 0
%...
%...   16A +   B +   C + 16D = 0
%...
%...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
%...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
%...                                                            4X
%...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
%...               X
%...                                                          4
%...  UI  = (1/12DX)(UI-2 - 8UI-1 + 0UI + 8UI+1 - UI+2) + O(DX )  (3)
%...    X
%...
%...     (4)  INTERIOR POINT, I = N-1
%...
%...                                              2               3
%...  A(UN-4 = UN-1 + UN-1  (-3DX)  + UN-1  (-3DX)  + UN-1  (-3DX)
%...                      X    1F         2X   2F         3X   3F
%...
%...                       4               5               6
%...         + UN-1  (-3DX)  + UN-1  (-3DX)  + UN-1  (-3DX)  + ...
%...               4X   4F         5X   5F         6X   6F
%...
%...                                              2               3
%...  B(UN-3 = UN-1 + UN-1  (-2DX)  + UN-1  (-2DX)  + UN-1  (-2DX)
%...                      X    1F         2X   2F         3X   3F
%...
%...                       4               5               6
%...         + UN-1  (-2DX)  + UN-1  (-2DX)  + UN-1  (-2DX)  + ...
%...               4X   4F         5X   5F         6X   6F
%...
%...                                              2               3
%...  C(UN-2 = UN-1 + UN-1  ( -DX)  + UN-1  (- -X)  + UN-1  ( -DX)
%...                      X    1F         2X   2F         3X   3F
%...
%...                       4               5               6
%...         + UN-1  ( -DX)  + UN-1  ( -DX)  + UN-1  ( -DX)  + ...
%...               4X   4F         5X   5F         6X   6F
%...
%...                                              2               3
%...  D(UN   = UN-1 + UN-1  (  DX)  + UN-1  (  DX)  + UN-1  (  DX)
%...                      X    1F         2X   2F         3X   3F
%...
%...                       4               5               6
%...         + UN-1  (  DX)  + UN-1  (  DX)  + UN-1  (  DX)  + ...
%...               4X   4F         5X   5F         6X   6F
%...
%...  -3A -  2B -   C +   D = 1
%...
%...   9A +  4B +   C +   D = 0
%...
%... -27A -  8B -   C +   D = 0
%...
%...  81A + 16B +   C +   D = 0
%...
%...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
%...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
%...                                                            4X
%...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
%...               X
%...                                                                4
%...  UN-1  = (1/12DX)(-UN-4 + 6UN-3 - 18UN-2 + 10UN-1 + 3UN) + O(DX )
%...      X
%...                                                              (4)
%...
%...    (5)  RIGHT END, POINT I = N
%...
%...                                       2             3
%...  A(UN-4 = UN + UN (-4DX)  + UN  (-4DX)  + UN  (-4DX)
%...                  X   1F       2X   2F       3X   3F
%...
%...                         4             5             6
%...             + UN  (-4DX)  + UN  (-4DX)  + UN  (-4DX)  + ...)
%...                 4X   4F       5X   5F       6X   6F
%...
%...                                       2             3
%...  B(UN-3 = UN + UN (-3DX)  + UN  (-3DX)  + UN  (-3DX)
%...                  X   1F       2X   2F       3X   3F
%...
%...                         4             5             6
%...             + UN  (-3DX)  + UN  (-3DX)  + UN  (-3DX)  + ...)
%...                 4X   4F       5X   5F       6X   6F
%...
%...                                       2             3
%...  C(UN-2 = UN + UN (-2DX)  + UN  (-2DX)  + UN  (-2DX)
%...                  X   1F       2X   2F       3X   3F
%...
%...                         4             5             6
%...             + UN  (-2DX)  + UN  (-2DX)  + UN  (-2DX)  + ...)
%...                 4X   4F       5X   5F       6X   6F
%...
%...                                       2             3
%...  D(UN-1 = UN + UN ( -DX)  + UN  ( -DX)  + UN  ( -DX)
%...                  X   1F       2X   2F       3X   3F
%...
%...                         4             5             6
%...             + UN  ( -DX)  + UN  ( -DX)  + UN  ( -DX)  + ...)
%...                 4X   4F       5X   5F       6X   6F
%...
%...   -4A -  3B -  2C -   D = 1
%...
%...   16A +  9B +  4C +   D = 0
%...
%...  -64A - 27B -  8C -   D = 0
%...
%...  256A + 81B + 16C +   D = 0
%...
%...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
%...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
%...                                                            4X
%...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
%...               X
%...                                                                4
%...  UN  = (1/12DX)(3UN-4 - 16UN-3 + 36UN-2 - 48UN-1 + 25UN) + O(DX )
%...    X
%...                                                              (5)
%...
%...  THE WEIGHTING COEFFICIENTS FOR EQUATIONS (1) TO (5) CAN BE
%...  SUMMARIZED AS
%...
%...             -25   48  -36   16   -3
%...
%...              -3  -10   18   -6    1
%...
%...       1/12    1   -8    0    8   -1
%...
%...              -1    6  -18   10    3
%...
%...               3  -16   36  -48   25
%...
%...  WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 4, M =
%...  1, P = 0, 1, 2, 3, 4 (BICKLEY, W. G., FORMULAE FOR NUMERICAL
%...  DIFFERENTIATION, MATH. GAZ., VOL. 25, 1941.  NOTE - THE BICKLEY
%...  COEFFICIENTS HAVE BEEN DIVIDED BY A COMMON FACTOR OF TWO).
%...
%...  EQUATIONS (1) TO (5) CAN NOW BE PROGRAMMED TO GENERATE THE
%...  DERIVATIVE U (X) OF FUNCTION u(X).
%...              X
%...
%...  COMPUTE THE SPATIAL INCREMENT
      dx=(xu-xl)/(n-1);
      r4fdx=1./(12.*dx);
      nm2=n-2;
%...
%...  EQUATION (1) (NOTE - THE RHS OF EQUATIONS (1), (2), (3), (4)
%...  AND (5) HAVE BEEN FORMATTED SO THAT THE NUMERICAL WEIGHTING
%...  COEFFICIENTS CAN BE MORE EASILY ASSOCIATED WITH THE BICKLEY
%...  MATRIX ABOVE)
      ux(  1)=r4fdx*...
      ( -25.*u(  1) +48.*u(  2) -36.*u(  3) +16.*u(  4)  -3.*u(  5));
%...
%...  EQUATION (2)
      ux(  2)=r4fdx*...
      (  -3.*u(  1) -10.*u(  2) +18.*u(  3)  -6.*u(  4)  +1.*u(  5));
%...
%...  EQUATION (3)
      for i=3:nm2
      ux(  i)=r4fdx*...
      (  +1.*u(i-2)  -8.*u(i-1)  +0.*u(  i)  +8.*u(i+1)  -1.*u(i+2));
      end
%...
%...  EQUATION (4)
      ux(n-1)=r4fdx*...
      (  -1.*u(n-4)  +6.*u(n-3) -18.*u(n-2) +10.*u(n-1)  +3.*u(  n));
%...
%...  EQUATION (5)
      ux(  n)=r4fdx*...
      (   3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u(  n));
