%DECK DSS044
      function [uxx]=dss044(xl,xu,n,u,ux,nl,nu)
%...
%...  FUNCTION DSS044 COMPUTES A FOURTH-ORDER APPROXIMATION OF A
%...  SECOND-ORDER DERIVATIVE, WITH OR WITHOUT THE NORMAL DERIVATIVE
%...  AT THE BOUNDARY.
%...
%...  ARGUMENT LIST
%...
%...     XL      LEFT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT)
%...
%...     XU      RIGHT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT)
%...
%...     N       NUMBER OF SPATIAL GRID POINTS, INCLUDING THE END
%...             POINTS (INPUT)
%...
%...     U       ONE-DIMENSIONAL ARRAY OF THE DEPENDENT VARIABLE TO BE
%...             DIFFERENTIATED (INPUT)
%...
%...     UX      ONE-DIMENSIONAL ARRAY OF THE FIRST DERIVATIVE OF U.
%...             THE END VALUES OF UX, UX(1) AND UX(N), ARE USED IN
%...             NEUMANN BOUNDARY CONDITIONS AT X = XL AND X = XU,
%...             DEPENDING ON THE ARGUMENTS NL AND NU (SEE THE DE-
%...             SCRIPTION OF NL AND NU BELOW)
%...
%...     UXX     ONE-DIMENSIONAL ARRAY OF THE SECOND DERIVATIVE OF U
%...             (OUTPUT)
%...
%...     NL      INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT
%...             X = XL (INPUT).  THE ALLOWABLE VALUES ARE
%...
%...                1 - DIRICHLET BOUNDARY CONDITION AT X = XL
%...                    (UX(1) IS NOT USED)
%...
%...                2 - NEUMANN BOUNDARY CONDITION AT X = XL
%...                    (UX(1) IS USED)
%...
%...     NU      INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT
%...             X = XU (INPUT).  THE ALLOWABLE VALUES ARE
%...
%...                1 - DIRICHLET BOUNDARY CONDITION AT X = XU
%...                    (UX(N) IS NOT USED)
%...
%...                2 - NEUMANN BOUNDARY CONDITION AT X = XU
%...                    (UX(N) IS USED)
%...
%...  THE FOLLOWING DERIVATION WAS COMPLETED BY W. E. SCHIESSER, DEPTS
%...  OF CHE AND CSEE, LEHIGH UNIVERSITY, BETHLEHEM, PA 18015, USA, ON
%...  DECEMBER 15, 1986.  ADDITIONAL DETAILS ARE GIVEN IN FUNCTION
%...  DSS042.
%...
%...  ******************************************************************
%...
%...  (1)  UXX AT THE INTERIOR POINTS 3, 4,..., N-2
%...
%...  TO DEVELOP A SET OF FOURTH-ORDER CORRECT DIFFERENTIATION FORMULAS
%...  FOR THE SECOND DERIVATIVE UXX, WE CONSIDER FIRST THE INTERIOR
%...  GRID POINTS AT WHICH A SYMMETRIC FORMULA CAN BE USED.
%...
%...  IF WE CONSIDER A FORMULA OF THE FORM
%...
%...     A*U(I-2) + B*U(I-1) + E*U(I) + C*U(I+1) + D*U(I+2)
%...
%...  TAYLOR SERIES EXPANSIONS OF U(I-2), U(I-1), U(I+1) AND U(I+2)
%...  CAN BE SUBSTITUTED INTO THIS FORMULA.  WE THEN CONSIDER THE
%...  LINEAR ALBEGRAIC EQUATIONS RELATING A, B, C AND D WHICH WILL
%...  RETAIN CERTAIN TERMS, I.E., UXX, AND DROP OTHERS, E.G., UXXX,
%...  UXXXX AND UXXXXX.
%...
%...  THUS, FOR GRID POINTS 3, 4,..., N-2
%...
%...     TO RETAIN UXX
%...
%...        4*A + B + C +  4*D = 2                              (1)
%...
%...     TO DROP UXXX
%...
%...       -8*A - B + C +  8*D = 0                              (2)
%...
%...     TO DROP UXXXX
%...
%...       16*A + B + C + 16*D = 0                              (3)
%...
%...     TO DROP UXXXXX
%...
%...      -32*A - B + C + 32*D = 0                              (4)
%...
%...  EQUATIONS (1) TO (4) CAN BE SOLVED FOR A, B, C AND D.  IF EQUA-
%...  TION (1) IS ADDED TO EQUATION (2)
%...
%...        -4*A + 2*C + 12*D = 2                               (5)
%...
%...  IF EQUATION (1) IS SUBTRACTED FROM EQUATION (3)
%...
%...        12*A + 12*D = -2                                    (6)
%...
%...  IF EQUATION (1) IS ADDED TO EQUATION (4)
%...
%...        -28*A + 2*C + 36*D = 2                              (7)
%...
%...  EQUATIONS (5) TO (7) CAN BE SOLVED FOR A, C AND D.  IF EQUATION
%...  (5) IS SUBTRACTED FROM EQUATION (7), AND THE RESULT COMBINED
%...  WITH EQUATION (6)
%...
%...         12*A + 12*D = -2                                   (6)
%...
%...        -24*A + 24*D = 0                                    (8)
%...
%...  EQUATIONS (6) AND (8) CAN BE SOLVED FOR A AND D.  FROM (8), A
%...  = D.  FROM EQUATION (6), A = -1/12 AND D = -1/12.  THEN, FROM
%...  EQUATION (5), C = 4/3, AND FROM EQUATION (1), B = 4/3.
%...
%...  THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS
%...
%...     (-1/12)*U(I-2) +   (4/3)*U(I-1) +
%...
%...       (4/3)*U(I+1) + (-1/12)*U(I+2)
%...
%...     (-1/12 + 4/3 - 1/12 + 4/3)*U(I) + UXX(I)*(DX**2) + O(DX**6)
%...
%...  OR
%...
%...     UXX(I) = (1/(12*DX**2))*(-1*U(I-2) + 16*U(I-1)
%...
%...                  - 30*U(I) + 16*U(I+1) -  1*U(I+2)         (9)
%...
%...                  + O(DX**4)
%...
%...  NOTE THAT THE UX TERM DROPS OUT, I.E., THE BASIC EQUATION IS
%...
%...        -2*A - B + C + 2*D =
%...
%...        -2*(-1/12) - (4/3) + (4/3) + 2*(-1/12) = 0
%...
%...  EQUATION (9) WAS OBTAINED BY DROPPING ALL TERMS IN THE UNDERLYING
%...  TAYLOR SERIES UP TO AND INCLUDING THE FIFTH DERIVATIVE, UXXXXX.
%...  THUS, EQUATION (9) IS EXACT FOR POLYNOMIALS UP TO AND INCLUDING
%...  FIFTH ORDER.  THIS CAN BE CHECKED BY SUBSTITUTING THE FUNCTIONS
%...  1, X, X**2, X**3, X**4 AND X**5 IN EQUATION (9) AND COMPUTING THE
%...  CORRESPONDING DERIVATIVES FOR COMPARISON WITH THE KNOWN SECOND
%...  DERIVATIVES.  THIS IS DONE FOR 1 MERELY BY SUMMING THE WEIGHTING
%...  COEFFICIENTS IN EQUATION (9), WHICH SHOULD SUM TO ZERO, I.E.,
%...  -1 + 16 - 30 + 16 -1 = 0.
%...
%...  FOR THE REMAINING FUNCTIONS, THE ALGEBRA IS RATHER INVOLVED, BUT
%...  THESE FUNCTIONS CAN BE CHECKED NUMERICALLY, I.E., NUMERICAL VALUES
%...  OF X**2, X**3, X**4 AND X**5 CAN BE SUBSTITUTED IN EQUATION (9)
%...  AND THE COMPUTED DERIVATIVES CAN BE COMPARED WITH THE KNOW NUMERI-
%...  CAL SECOND DERIVATIVES.  THIS IS NOT A PROOF OF CORRECTNESS OF
%...  EQUATION (9), BUT WOULD LIKELY DETECT ANY ERRORS IN EQUATION (9).
%...
%...  CHECK THE VALUES OF A, B, C AND D IN EQUATIONS (1) TO (4)
%...  CALL CHECK(1)
%...
%...  ******************************************************************
%...
%...  (2)  UXX AT THE INTERIOR POINTS I = 2 AND N-1
%...
%...  FOR GRID POINT 2, WE CONSIDER A FORMULA OF THE FORM
%...
%...     A*U(I-1) + F*U(I) + B*U(I+1) + C*U(I+2) + D*U(I+3) + E*U(I+4)
%...
%...  TAYLOR SERIES EXPANSIONS OF U(I-1), U(I+1), U(I+2), U(I+3) AND
%...  U(I+4) WHEN SUBSTITUTED INTO THIS FORMULA GIVE LINEAR ALGEBRAIC
%...  EQUATIONS RELATING A, B, C, D AND E.
%...
%...     TO DROP UX
%...
%...        -A + B + 2*C + 3*D + 4*E = 0                       (10)
%...
%...     TO RETAIN UXX
%...
%...        A + B + 4*C + 9*D + 16*E = 2                       (11)
%...
%...     TO DROP UXXX
%...
%...       -A + B + 8*C + 27*D + 64*E = 0                      (12)
%...
%...     TO DROP UXXXX
%...
%...        A + B + 16*C + 81*D + 256*E = 0                    (13)
%...
%...     TO DROP UXXXXX
%...
%...       -A + B + 32*C + 243*D + 1024*E = 0                  (14)
%...
%...  EQUATIONS (11), (12), (13) AND (14) CAN BE SOLVED FOR A, B, C,
%...  D AND E.  IF EQUATION (10) IS ADDED TO EQUATION (11)
%...
%...     2*B + 6*C + 12*D +20*E = 2                            (15)
%...
%...  IF EQUATION (10) IS SUBTRACTED FROM EQUATION (12)
%...
%...     6*C + 24*D + 60*E = 0                                 (16)
%...
%...  IF EQUATION (10) IS ADDED TO EQUATION (13)
%...
%...     2*B + 18*C + 84*D + 260*E = 0                         (17)
%...
%...  IF EQUATION (10) IS SUBTRACTED FROM EQUATION (14)
%...
%...     30*C + 240*D + 1020*E = 0                             (18)
%...
%...  EQUATIONS (15), (16), (17) AND (18) CAN BE SOLVED FOR B, C, D
%...  AND E.
%...
%...     6*C + 24*D + 60*E = 0                                 (16)
%...
%...  IF EQUATION (15) IS SUBTRACTED FROM EQUATION (17)
%...
%...     12*C + 72*D + 240*E = -2                              (19)
%...
%...     30*C + 240*D + 1020*E = 0                             (18)
%...
%...  EQUATIONS (16), (18) AND (19) CAN BE SOLVED FOR C, D AND E.  IF
%...  TWO TIMES EQUATION (16) IS SUBTRACTED FROM EQUATION (19),
%...
%...     24*D + 120*E = -2                                     (20)
%...
%...  IF FIVE TIMES EQUATION (16) IS SUBTRACTED FROM EQUATION (18),
%...
%...     120*D + 720*E = 0                                     (21)
%...
%...  EQUATIONS (20) AND (21) CAN BE SOLVED FOR D AND E.  FROM (21),
%...  E = (-1/6)*D.  SUBSTITUTION IN EQUATION (20) GIVES D = -1/2.
%...  THUS, E = 1/12.  FROM EQUATION (16), C = 7/6.  FROM EQUATION
%...  (15), B = -1/3.  FROM EQUATION (10), A = 5/6.
%...
%...  THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS
%...
%...  (5/6)*U(I-1) + (-1/3)*U(I+1) + (7/6)*U(I+2) + (-1/2)*U(I+3)
%...
%...  + (1/12)*U(I+4) = (5/6 - 1/3 + 7/6 - 1/2 + 1/12)*U(I)
%...
%...  + UXX*(DX**2) + O(DX**6)
%...
%...  OR
%...
%...  UXX(I) = (1/12*DX**2)*(10*U(I-1) - 15*U(I) - 4*U(I+1)
%...                                                           (22)
%...         + 14*U(I+2) - 6*U(I+3) + 1*U(I+4)) + O(DX**4)
%...
%...  EQUATION (22) WILL BE APPLIED AT I = 2 AND N-1.  THUS
%...
%...  UXX(2) = (1/12*DX**2)*(10*U(1) - 15*U(2) - 4*U(3)
%...                                                           (23)
%...         + 14*U(4) - 6*U(5) + 1*U(6)) + O(DX**4)
%...
%...  UXX(N-1) = (1/12*DX**2)*(10*U(N) - 15*U(N-1) - 4*U(N-2)
%...                                                           (24)
%...           + 14*U(N-3) - 6*U(N-4) + 1*U(N-5)) + O(DX**4)
%...
%...  CHECK THE VALUES OF A, B, C, D AND E IN EQUATIONS (10) TO (14)
%...  CALL CHECK(2)
%...
%...  ******************************************************************
%...
%...  (3)  UXX AT THE BOUNDARY POINTS 1 AND N
%...
%...  FINALLY, FOR GRID POINT 1, AN APPROXIMATION WITH A NEUMANN BOUND-
%...  ARY CONDITON OF THE FORM
%...
%...     A*U(I+1) + B*U(I+2) + C*U(I+3) + D*U(I+4) + E*UX(I) + F*U(I)
%...
%...  WILL BE USED.  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
%...
%...     TO DROP UX
%...
%...        A + 2*B + 3*C + 4*D + E = 0                        (25)
%...
%...     TO RETAIN UXX
%...
%...        A + 4*B + 9*C + 16*D = 2                           (26)
%...
%...     TO DROP UXXX
%...
%...        A + 8*B + 27*C + 64*D = 0                          (27)
%...
%...     TO DROP UXXXX
%...
%...        A + 16*B + 81*C + 256*D = 0                        (28)
%...
%...     TO DROP UXXXXX
%...
%...        A + 32*B + 243*C + 1024*D = 0                      (29)
%...
%...  EQUATIONS (25) TO (29) CAN BE SOLVED FOR A, B, C, D AND E.  IF
%...
%...  EQUATION (26) IS SUBTRACTED FROM EQUATIONS (27), (28) AND (29),
%...
%...     4*B + 18*C + 48*D = -2                                (30)
%...
%...     12*B + 72*C + 240*D = -2                              (31)
%...
%...     28*B + 234*C + 1008*D = -2                            (32)
%...
%...  EQUATIONS (30), (31) AND (32) CAN BE SOLVED FOR B, C AND D
%...
%...     18*C + 96*D = 4                                       (33)
%...
%...     108*C + 672*D = 12                                    (34)
%...
%...  EQUATIONS (3) AND (34) CAN BE SOLVED FOR C AND D, C = 8/9, D =
%...  -1/8.
%...
%...  FROM EQUATION (30), B = -3.  FROM EQUATION (26), A = 8.  FROM
%...  EQUATION (25), E = -25/6.
%...
%...  THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS
%...
%...  8*U(I+1) - 3*U(I+2) + (8/9)*U(I+3) - (1/8)*U(I+4)
%...
%...  - (25/6)*UX(I)*DX
%...
%...  = (8 - 3 + (8/9) - (1/8))*U(I) + UXX*(DX**2) + O(DX**6)
%...
%...  OR
%...
%...  UXX(I) = (1/12*DX**2)*((-415/6)*U(I) + 96*U(I+1) - 36*U(I+2)
%...                                                                (35)
%...  + (32/3)*U(I+3) - (3/2)*U(I+4) - 50*UX(I)*DX) + O(DX**4)
%...
%...  EQUATION (35) WILL BE APPLIED AT I = 1 AND I = N
%...
%...  UXX(1) = (1/12*DX**2)*((-415/6)*U(1) + 96*U(2) - 36*U(3)
%...                                                                (36)
%...  + (32/3)*U(4) - (3/2)*U(5) - 50*UX(1)*DX) + O(DX**4)
%...
%...  UXX(N) = (1/12*DX**2)*((-415/6)*U(N) + 96*U(N-1) - 36*U(N-2)
%...                                                                (37)
%...  + (32/3)*U(N-3) - (3/2)*U(N-4) + 50*UX(N)*DX) + O(DX**4)
%...
%...  CHECK THE VALUES OF A, B, C AND D IN EQUATIONS (25) TO (29)
%...  CALL CHECK(3)
%...
%...  ALTERNATIVELY, FOR GRID POINT 1, AN APPROXIMATION WITH A DIRICHLET
%...  BOUNDARY CONDITION OF THE FORM
%...
%...  A*U(I+1) + B*U(I+2) + C*U(I+3) + D*U(I+4) + E*U(I+5) + F*U(I)
%...
%...  CAN BE USED.  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
%...
%...     TO DROP UX
%...
%...        A + 2*B + 3*C + 4*D + 5*E = 0                      (38)
%...
%...     TO RETAIN UXX
%...
%...        A + 4*B + 9*C + 16*D + 25*E = 2                    (39)
%...
%...     TO DROP UXXX
%...
%...        A + 8*B + 27*C + 64*D + 125*E = 0                  (40)
%...
%...     TO DROP UXXXX
%...
%...        A + 16*B + 81*C + 256*D + 625*E = 0                (41)
%...
%...     TO DROP UXXXXX
%...
%...        A + 32*B + 243*C + 1024*D + 3125*E = 0             (42)
%...
%...  EQUATIONS (38), (39), (40), (41) AMD (42) CAN BE SOLVED FOR A,
%...  B, C, D AND E.
%...
%...        2*B + 6*C + 12*D + 20*E = 2                        (43)
%...
%...        6*B + 24*C + 60*D + 120*E = 0                      (44)
%...
%...        14*B + 78*C + 252*D + 620*E = 0                    (45)
%...
%...        30*B + 240*C + 1020*D + 3120*E = 0                 (46)
%...
%...  EQUATIONS (43), (44), (45) AND (46) CAN BE SOLVED FOR B, C, D
%...  AND E
%...
%...        6*C + 24*D + 60*E = -6                             (47)
%...
%...        36*C + 168*D + 480*E = -14                         (48)
%...
%...        150*C + 840*D + 2820*E = -30                       (49)
%...
%...  EQUATIONS (47), (48) AND (49) CAN BE SOLVED FOR C, D AND E
%...
%...        24*D + 120*E = 22                                  (50)
%...
%...        240*D + 1320*E = 120                               (51)
%...
%...  FROM EQUATIONS (50) AND (51), D = 61/12, E = -5/6.  FROM EQUATION
%...  (47), C = -13.  FROM EQUATION (43), B = 107/6.  FROM EQUATION
%...  (38), A = -77/6.
%...
%...  THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS
%...
%...  (-77/6)*U(I+1) + (107/6)*U(I+2) - 13*U(I+3) + (61/12)*U(I+4)
%...
%...  - (5/6)*U(I+5) = (-77/6 + 107/6 - 13 + 61/12 - 5/6)*U(I) +
%...
%...  UXX(I)*(DX**2) + O(DX**6)
%...
%...  OR
%...
%...  UXX(I) = (1/12*DX**2)*(45*U(I) - 154*U(I+1) + 214*U(I+2)
%...                                                                (52)
%...         - 156*U(I+3) + 61*U(I+4) - 10*U(I+5)) + O(DX**4)
%...
%...  EQUATION (52) WILL BE APPLIED AT I = 1 AND I = N
%...
%...  UXX(1) = (1/12*DX**2)*(45*U(1) - 154*U(2) + 214*U(3)
%...                                                                (53)
%...         - 156*U(4) + 61*U(5) - 10*U(6)) + O(DX**4)
%...
%...  UXX(N) = (1/12*DX**2)*(45*U(N) - 154*U(N-1) + 214*U(N-2)
%...                                                                (54)
%...         -156*U(N-3) + 61*U(N-4) - 10*U(N-5)) + O(DX**4)
%...
%...  CHECK THE VALUES OF A, B, C, D, AND E IN EQUATIONS (38) TO (42)
%...  CALL CHECK(4)
%...
%...  ******************************************************************
%...
%...  GRID SPACING
      dx=(xu-xl)/(n-1);
%...
%...  1/(12*DX**2) FOR SUBSEQUENT USE
      r12dxs=1./(12.0*dx^2);
%...
%...  UXX AT THE LEFT BOUNDARY
%...
%...     WITHOUT UX (EQUATION (53))
         if nl==1
         uxx(1)=r12dxs*...
                       (   45.0*u(1)...
                         -154.0*u(2)...
                         +214.0*u(3)...
                         -156.0*u(4)...
                          +61.0*u(5)...
                          -10.0*u(6));
%...
%...     WITH UX (EQUATION (36))
         elseif nl==2
         uxx(1)=r12dxs*...
                       (-415.0e0/6.0*u(1)...
                               +96.0*u(2)...
                               -36.0*u(3)...
                         +32.0e0/3.0*u(4)...
                          -3.0e0/2.0*u(5)...
                               -50.0*ux(1)*dx);
         end
%...
%...  UXX AT THE RIGHT BOUNDARY
%...
%...     WITHOUT UX (EQUATION (54))
         if nu==1
         uxx(n)=r12dxs*...
                       (   45.0*u(n  )...
                         -154.0*u(n-1)...
                         +214.0*u(n-2)...
                         -156.0*u(n-3)...
                          +61.0*u(n-4)...
                          -10.0*u(n-5));
%...
%...     WITH UX (EQUATION (37))
         elseif nu==2
         uxx(n)=r12dxs*...
                       (-415.0e0/6.0*u(n  )...
                               +96.0*u(n-1)...
                               -36.0*u(n-2)...
                         +32.0e0/3.0*u(n-3)...
                          -3.0e0/2.0*u(n-4)...
                               +50.0*ux(n  )*dx);
         end
%...
%...  UXX AT THE INTERIOR GRID POINTS
%...
%...     I = 2 (EQUATION (23))
         uxx(2)=r12dxs*...
                       (   10.0*u(1)...
                          -15.0*u(2)...
                           -4.0*u(3)...
                          +14.0*u(4)...
                           -6.0*u(5)...
                           +1.0*u(6));
%...
%...     I = N-1 (EQUATION (24))
         uxx(n-1)=r12dxs*...
                         ( 10.0*u(n  )...
                          -15.0*u(n-1)...
                           -4.0*u(n-2)...
                          +14.0*u(n-3)...
                           -6.0*u(n-4)...
                           +1.0*u(n-5));
%...
%...     I = 3, 4,..., N-2 (EQUATION (9))
         for i=3:n-2
         uxx(i)=r12dxs*...
                       (   -1.0*u(i-2)...
                          +16.0*u(i-1)...
                          -30.0*u(i  )...
                          +16.0*u(i+1)...
                           -1.0*u(i+2));
         end
