%DECK DSS002
      function [ux]=dss002(xl,xu,n,u)
%...
%...  FUNCTION DSS002 COMPUTES THE FIRST DERIVATIVE, U , OF A
%...                                                  X
%...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU
%...
%...  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)
%...
%...  FUNCTION DSS002 COMPUTES THE FIRST DERIVATIVE, U , OF A
%...                                                  X
%...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM THE
%...  CLASSICAL THREE-POINT, SECOND-ORDER FINITE DIFFERENCE APPROXI-
%...  TIONS
%...
%...                                       2
%...  U1  = (1/2DX)(-3U1 + 4U2 - U3) + O(DX ) (LEFT BOUNDARY,     (1)
%...    X                                         X = XL)
%...
%...                                   2
%...  UI  = (1/2DX)(UI+1 - UI-1) + O(DX ) (INTERIOR POINT,        (2)
%...    X                                   X NE XL, XU)
%...
%...                                          2
%...  UN  = (1/2DX)(3UN - 4UN-1 + UN-2) + O(DX ) (RIGHT BOUNDARY, (3)
%...    X                                            X = XU)
%...
%...  EQUATIONS (1) TO (3) APPLY OVER A GRID IN X WITH CORRESPONDING
%...  VALUES OF THE FUNCTION U(X) REPRESENTED AS
%...
%...   U1      U2       U3         UI        UN-2      UN-1    UN
%...
%...  X=XL  X=XL+DX  X=XL+2DX ... X=XI ... X=XU-2DX  X=XU-DX  X=XU
%...
%...  THE ORIGIN OF EQUATIONS (1) TO (3) IS OUTLINED BELOW.
%...
%...  CONSIDER THE FOLLOWING POLYNOMIAL IN X OF ARBITRARY ORDER
%...
%...                                     2             3
%...  U(X) = A0 + A1(X - X0) + A2(X - X0)  + A3(X - X0)  + ....   (4)
%...
%...  WE SEEK THE VALUES OF THE COEFFICIENTS A0, A1, A2, ... FOR A
%...  PARTICULAR FUNCTION U(X).  IF X = X0 IS SUBSTITUTED IN EQUATION
%...  (4), WE HAVE IMMEDIATELY A0 = U(X0).  NEXT, IF EQUATION (4) IS
%...  DIFFERENTIATED WITH RESPECT TO X,
%...
%...                                                   2
%...  DU(X)/DX = U (X) = A1 + 2A2(X - X0) + 3A3(X - X0)  + ...    (5)
%...              X
%...
%...  AGAIN, WITH X = X0, A1 = DU(X0)/DX = U (X0).  DIFFERENTIATION
%...                                        X
%...  OF EQUATION (5) IN TURN GIVES
%...
%...  D2U(X)/DX2 = U  (X) = 2A2 + 6A3(X - X0) + ...
%...                2X
%...
%...  AND FOR X = X0, A2 = U  (X0)/2F (2F = 1*2, I.E., 2 FACTORIAL).
%...                        2X
%...
%...  WE CAN CONTINUE THIS PROCESS OF DIFFERENTIATION FOLLOWED BY THE
%...  SUBSTITUTION X = X0 TO OBTAIN THE SUCCESSIVE COEFFICIENTS IN
%...  EQUATION (4), A3, A4, ...  FINALLY, SUBSTITUTION OF THESE CO-
%...  EFFICIENTS IN EQUATION (4) GIVES
%...
%...                                                 2
%...  U(X) = U(X0) + U (X0)(X - X0) + U  (X0)(X - X0)  +
%...                  X       1F       2X       2F
%...                                                              (6)
%...                                3                  4
%...                 U  (X0)(X - X0)  + U  (X0)(X - X0)  + ...
%...                  3X       3F        4X       4F
%...
%...  THE CORRESPONDENCE BETWEEN EQUATION (6) AND THE WELL-KNOWN
%...  TAYLOR SERIES SHOULD BE CLEAR.  THUS THE EXPANSION OF A
%...  FUNCTION, U(X), AROUND A NEIGHBORING POINT X0 IN TERMS OF U(X0)
%...  AND THE DERIVATIVES OF U(X) AT X = X0 IS EQUIVALENT TO APPROXI-
%...  MATING U(X) NEAR X0 BY A POLYNOMIAL.
%...
%...  EQUATION (6) IS THE STARTING POINT FOR THE DERIVATION OF THE
%...  CLASSICAL FINITE DIFFERENCE APPROXIMATIONS OF DERIVATIVES SUCH
%...  AS THE THREE-POINT FORMULAS OF EQUATIONS (1), (2) AND (3).  WE
%...  WILL NOW CONSIDER THE DERIVATION OF THESE THREE-POINT FORMULAS
%...  IN A STANDARD FORMAT WHICH CAN THEN BE EXTENDED TO HIGHER
%...  MULTI-POINT FORMULAS IN OTHER SUBROUTINES, E.G., FIVE-POINT
%...  FORMULAS IN SUBROUTINE DSS004.
%...
%...  THREE-POINT FORMULAS
%...
%...     (1)  LEFT END, POINT I = 1
%...
%...  IF EQUATION (6) IS WRITTEN AROUND THE POINT X = XL FOR X = XL +
%...  DX AND X = XL + 2DX, FOR WHICH THE CORRESPONDING VALUES OF U(X)
%...  ARE U1, U2 AND U3 (U1 AND U2 ARE SEPARATED WITH RESPECT TO X BY
%...  DISTANCE DX AS ARE U2 AND U3, I.E., WE ASSUME A UNIFORM GRID
%...  SPACING, DX, FOR INDEPENDENT VARIABLE X)
%...
%...                                2            3
%...  U2 = U1 + U1 ( DX) + U1  ( DX)  + U1  ( DX)  + ...          (7)
%...              X  1F      2X  2F       3X  3F
%...
%...                                2            3
%...  U3 = U1 + U1 (2DX) + U1  (2DX)  + U1  (2DX)  + ...          (8)
%...              X  1F      2X  2F       3X  3F
%...
%...  WE CAN NOW TAKE A LINEAR COMBINATION OF EQUATIONS (7) AND (8)
%...  BY FIRST MULTIPLYING EQUATION (7) BY A CONSTANT, A, AND EQUA-
%...  TION (8) BY CONSTANT B
%...
%...                                  2           3
%...  A(U2 = U1 + U1 ( DX) + U1  ( DX) + U1  ( DX) + ...)         (9)
%...                X  1F      2X  2F      3X  3F
%...
%...                                  2           3
%...  B(U3 = U1 + U1 (2DX) + U1  (2DX) + U1  (2DX) + ...)        (10)
%...                X  1F      2X  2F      3X  3F
%...
%...  CONSTANTS A AND B ARE THEN SELECTED SO THAT THE COEFFICIENTS OF
%...  THE U1  TERMS SUM TO ONE (SINCE WE ARE INTERESTED IN OBTAINING
%...        X
%...  A FINITE DIFFERENCE APPROXIMATION FOR THIS FIRST DERIVATIVE).
%...  ALSO, WE SELECT A AND B SO THAT THE COEFFICIENTS OF THE U1
%...                                                            2X
%...  TERMS SUM TO ZERO IN ORDER TO DROP OUT THE CONTRIBUTION OF THIS
%...  SECOND DERIVATIVE (THE BASIC IDEA IS TO DROP OUT AS MANY OF THE
%...  DERIVATIVES AS POSSIBLE IN THE TAYLOR SERIES BEYOND THE DERI-
%...  VATIVE OF INTEREST, IN THIS CASE U1 , IN ORDER TO PRODUCE A
%...                                     X
%...  FINITE DIFFERENCE APPROXIMATION FOR THE DERIVATIVE OF MAXIMUM
%...  ACCURACY).  IN THIS CASE WE HAVE ONLY TWO CONSTANTS, A AND B,
%...  TO SELECT SO WE CAN DROP OUT ONLY THE SECOND DERIVATIVE, U1  ,
%...                                                             2X
%...  IN THE TAYLOR SERIES (IN ADDITION TO RETAINING THE FIRST DERI-
%...  VATIVE).  THIS PROCEDURE LEADS TO TWO LINEAR ALGEBRAIC EQUA-
%...  TIONS IN THE TWO CONSTANTS
%...
%...  A + 2B = 1
%...
%...  A + 4B = 0
%...
%...  SOLUTION OF THESE EQUATIONS FOR A AND B GIVES
%...
%...  A = 2, B = -1/2
%...
%...  SOLUTION OF EQUATIONS (9) AND (10) FOR U1  WITH THESE VALUES OF
%...  A AND B GIVES EQUATION (1)               X
%...
%...                                      2
%...  U1 = (1/2DX)(-3U1 + 4U2 - U3) + O(DX )                      (1)
%...    X
%...               2
%...  THE TERM O(DX ) INDICATES A PRINCIPAL ERROR TERM DUE TO TRUNCA-
%...                                                2
%...  TION OF THE TAYLOR SERIES WHICH IS OF ORDER DX .  THIS TERM IN
%...                    2
%...  FACT EQUALS U1  DX /3F, WHICH IS EASILY OBTAINED IN DERIVING
%...                3X
%...  EQUATION (1).
%...
%...  THIS SAME BASIC PROCEDURE CAN NOW BE APPLIED TO THE DERIVATION
%...  OF EQUATIONS (2) AND (3).
%...
%...     (2)  INTERIOR POINT I
%...
%...                                    2           3
%...  A(UI-1 = UI + UI (-DX) + UI  (-DX) + UI  (-DX) + ...)
%...                  X  1F      2X  2F      3X  3F
%...
%...                                    2           3
%...  B(UI+1 = UI + UI ( DX) + UI  ( DX) + UI  ( DX) + ...)
%...                  X  1F      2X  2F      3X  3F
%...
%...  -A + B = 1
%...
%...   A + B = 0
%...
%...  A = 1/2, B = -1/2
%...                                   2
%...  UI  = (1/2DX)(UI+1 - UI-1) + O(DX )                         (2)
%...    X
%...
%...     (3)  RIGHT END, POINT I = N
%...
%...                                      2            3
%...  A(UN-2 = UN + UN (-2DX) + UN  (-2DX) + UN  (-2DX) + ...)
%...                  X   1F      2X   2F      3X   3F
%...
%...                                      2            3
%...  B(UN-1 = UN + UN ( -DX) + UN  ( -DX) + UN  ( -DX) + ...)
%...                  X   1F      2X   2F      3X   3F
%...
%...  -2A - B = 1
%...
%...   4A + B = 0
%...
%...   A = -2, B = 1/2
%...                                          2
%...  UN  = (1/2DX)(3UN - 4UN-1 + UN-2) + O(DX )                  (3)
%...    X
%...
%...  THE WEIGHTING COEFFICIENTS FOR EQUATIONS (1), (2) AND (3) CAN
%...  BE SUMMARIZED AS
%...
%...          -3   4  -1
%...
%...     1/2  -1   0   1
%...
%...           1  -4   3
%...
%...  WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 2, M =
%...  1, P = 0, 1, 2 (BICKLEY, W. G., FORMULAE FOR NUMERICAL DIFFER-
%...  ENTIATION, MATH. GAZ., VOL. 25, 1941).
%...
%...  EQUATIONS (1), (2) AND (3) CAN NOW BE PROGRAMMED TO GENERATE
%...  THE DERIVATIVE U (X) OF U(X).
%...                  X
%...
%...  COMPUTE THE SPATIAL INCREMENT
      dx=(xu-xl)/(n-1);
      r2fdx=1./(2.*dx);
      nm1=n-1;
%...
%...  EQUATION (1) (NOTE - THE RHS OF THE FINITE DIFFERENCE APPROXI-
%...  TIONS, EQUATIONS (1), (2) AND (3) HAVE BEEN FORMATTED SO THAT
%...  THE NUMERICAL WEIGHTING COEFFICIENTS CAN BE MORE EASILY ASSOCI-
%...  ATED WITH THE BICKLEY MATRIX LISTED ABOVE)
      ux(1)=r2fdx*...
      (     -3.   *u(  1)     +4.   *u(  2)     -1.   *u(  3));
%...
%...  EQUATION (2)
      for i=2:nm1
      ux(i)=r2fdx*...
      (     -1.   *u(i-1)     +0.   *u(  i)     +1.   *u(i+1));
      end
%...
%...  EQUATION (3)
      ux(n)=r2fdx*...
      (      1.   *u(n-2)     -4.   *u(n-1)     +3.   *u(  n));
