dss042

%  File: dss042.m
%
   function [uxx]=dss042(xl,xu,n,u,ux,nl,nu)
%
%  Function dss042 computes a second-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 is for a set of second-order, four-point
%  approximations for a second derivative that can be used at the
%  boundaries of a spatial domain.  These approximations have the
%  features
%
%     (1)  Only interior and boundary points are used (i.e., no
%          fictitious points are used)
%
%     (2)  The normal derivative at the boundary is included as part
%          of the approximation for the second derivative
%
%     (3)  Approximations for the boundary conditions are not used.
%
%  The derivation is by Professor Gilbert A. Stengle, Department of
%  Mathematics, Lehigh University, Bethlehem, pa 18015, and was done
%  on December 7, 1985.
%
%  For an approximation at the left boundary, involving the points
%  i, i+1, i+2 and i+3, consider the following Taylor series expan-
%  sions
%
%                  ux(i)( dx)   uxx(i)( dx)**2   uxxx(i)( dx)**3
%  u(i+1) = u(i) + ---------- + -------------- + --------------- +...
%                       1             2                 6
%
%
%                  ux(i)(2dx)   uxx(i)(2dx)**2   uxxx(i)(2dx)**3
%  u(i+2) = u(i) + ---------- + -------------- + --------------- +...
%                       1             2                 6
%
%  If we now form the following linear combination, involving con-
%  stants a, b, c and d to be determined, and use the preceding two
%  Taylor series,
%
%     a*u(i) + b*ux(i) + c*u(i+1) + d*u(i+2)
%
%  we have
%
%     a*u(i) + b*ux(i) + c*u(i+1) + d*u(i+2) =
%
%     (a + b + c + d)*u(i) +
%
%     (b + dx*c + 2*dx*d)*ux(i) +
%
%     (c*(dx**2)/2 + d*((2*dx)**2)/2)*uxx(i) +
%
%     (c*(dx**3)/6 + d*((2*dx)**3)/6)*uxxx(i) + O(dx**4)
%
%  The third derivative, uxxx(i), can be dropped by taking
%
%     c = -8*d
%
%  The second derivative, uxx(i), can be retained by taking
%
%     (dx**2)(c/2 + 2*d) = 1
%
%  which, when combined with the preceding result gives
%
%     d = -1/(2*(dx**2))
%
%     c = 4/(dx**2)
%
%  The first derivative, ux(i), can be dropped by taking
%
%     b + dx*c + 2*dx*d = 0
%
%  or
%
%     b = -dx*c - 2*dx*d = -4/dx - 2*dx*(-1/(2*(dx**2))) = -3/dx
%
%  Finally, u(i), can be dropped by taking
%
%     a = - c - d = 8*d - d = -7*d = -7/(2*(dx**2))
%
%  If we now solve for the derivative of interest, uxx(i),
%
%     uxx(i) = -7/(2(dx**2))*u(i) - 3/dx*ux(i)
%
%              + 8/(dx**2)*u(i+1) - 1/(2*(dx**2))u(i+2) + O(dx**2)
%
%       = (1/(2*(dx**2)))*(-u(i+2) + 8*u(i+1) - 7*u(i) - 6*dx*ux(i))
%
%         + O(dx**2)
%
%  which is the four-point, second-order approximation for the second
%  derivative, uxx(i), including the first derivative, ux(i).
%
%  Four checks of this approximation can easily be made for u(i) =
%  1, u(i) = x, u(i) = x**2 and u(i) = x**3
%
%     uxx(i) = (1/(2*(dx**2)))*(-1 + 8*1 - 7*1 - 6*dx*0) = 0
%
%     uxx(i) = (1/(2*(dx**2)))*(-(x + 2*dx) + 8*(x + dx)
%
%              -7*x - 6*dx*1) = 0
%
%     uxx(i) = (1/(2*(dx**2)))*(-(x + 2*dx)**2 + 8*(x + dx)**2
%
%            - 7*(x**2) - 6*dx*(2*x))
%
%             = (-  x**2 -  4*x*dx - 4*dx**2
%
%               + 8*x**2 + 16*x*dx + 8*dx**2
%
%               - 7*x**2 - 12*x*dx)/(2*(dx**2)) = 2
%
%     uxx(i) = (1/(2*(dx**2)))*(-(x + 2*dx)**3 + 8*(x + dx)**3
%
%            - 7*(x**3) - 6*dx*(3*x**2))
%
%            = (1/(2*(dx**2)))*(- x**3 - 6*dx*x**2 - 12*x*dx**2
%
%            - 8*dx**3 + 8*x**3 + 24*dx*x**2 + 24*x*dx**2 + 8*dx**3
%
%            - 7*x**3 - 18*dx*x**2)
%
%            = (1/(2*(dx**2)))*(12*x*dx**2) = 6*x
%
%  The preceding approximation for uxx(i) can be applied at the
%  left boundary value of x by taking i = 1.  An approximation at
%  the right boundary is obtained by taking dx = -dx and reversing
%  the subscripts in the preceding approximation, with i = n
%
%     uxx(i)
%
%       = (1/(2*(dx**2)))*(-u(i-2) + 8*u(i-1) - 7*u(i) + 6*dx*ux(i))
%
%         + O(dx**2)
%
%  To obtain approximations of the second derivative which do not
%  involve the first derivative, we take as the linear combination
%
%     a*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3) 
%
%  we have
%
%     a*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3) =
%
%     (a + b + c + d)*u(i)+
%
%     (dx*b + 2*dx*c + 4*dx*d)*ux(i)+
%
%     (b*(dx**2)/2 + c*((2*dx)**2)/2 + d*((3*dx)**2)/2)*uxx(i) +
%
%     (b*(dx**3)/6 + c*((2*dx)**3)/6 + d*((3*dx)**3)/6)*uxx(i) +
%
%     O(dx**4)
%
%  The third derivative, uxxx(i), can be dropped by taking
%
%     b + 8*c + 27*d = 0
%
%  The second derivative, uxx(i), can be retained by taking
%
%     (dx**2)*(b/2 + 2*c + (9/2)*d) = 1
%
%  The first derivative can be dropped by taking
%
%     b + 2*c + 3*d = 0
%
%  Solution of the preceding equations for c and d by elimination of
%  b gives
%
%     6*c + 24*d = 0
%
%     4*c + 18*d = -2/(dx**2)
%
%  Then, eliminating c, gives
%
%     (18 - 16)*d = -2/(dx**2)
%
%  or
%
%     d = -1/(dx**2)
%
%     c = (24/6)/(dx**2) = 4/(dx**2)
%
%     b = -8/(dx**2) + 3/(dx**2) = -5/(dx**2)
%
%  u(i) can be dropped by taking
%
%     a + b + c + d = 0
%
%  or
%
%     a = (5 - 4 + 1)/(dx**2) = 2/(dx**2)
%
%  If we now solve for the derivative of interest, uxx(i),
%
%     uxx(i) = (1/dx**2)*(2*u(i) - 5*u(i+1) + 4*u(i+2) - 1*u(i+3))
%
%            + O(dx**2)
%
%  Which is the four-point, second-order approximation for the second
%  derivative, uxx(i), without the first derivative, ux(i).
%
%  Four checks of this approximation can easily be made for u(i) =
%  1, u(i) = x, u(i) = x**2 and u(i) = x**3
%
%     uxx(i) = (1/dx**2)*(2 - 5 + 4 - 1) = 0
%
%     uxx(i) = (1/dx**2)*(2*x - 5*(x + dx) + 4*(x + 2*dx)
%
%              - 1*(x + 3*dx)) = 0
%
%     uxx(i) = (1/dx**2)*(2*x**2 - 5*(x + dx)**2 + 4*(x + 2*dx)**2
%
%              - 1*(x + 3*dx)**2) = 2
%
%     uxx(i) = (1/dx**2)*(2*x**3 - 5*(x + dx)**3 + 4*(x + 2*dx)**3
%
%            - 1*(x + 3*dx)**3)
%
%             = (1/dx**2)*(2*x**3 - 5*x**3 - 15*x*dx**2
%
%             - 15*dx*x**2 - 5*dx**3 + 4*x**3 + 24*dx*x**2
%
%             + 48*x*dx**2 + 32*dx**3 - x**3 - 9*dx*x**2
%
%             - 27*x*dx**2 - 27dx**3)
%
%             = (1/dx**2)*(6*x*dx**2) = 6*x
%
%  The preceding approximation for uxx(i) can be applied at the
%  left boundary value of x by taking i = 1.  An approximation at
%  the right boundary is obtained by taking dx = -dx and reversing
%  the subscripts in the preceding approximation, with i = n
%
%     uxx(i) = (1/dx**2)*(2*u(i) - 5*u(i-1) + 4*u(i-2) - 1*u(i-3))
%
%            + O(dx**2)
%
%  Grid spacing
   dx=(xu-xl)/(n-1);
%
%  uxx at the left boundary, without ux
   if nl==1
     uxx(1)=(( 2.)*u(  1)...
            +(-5.)*u(  2)...
            +( 4.)*u(  3)...
            +(-1.)*u(  4))/(dx^2);
%
%  uxx at the left boundary, including ux
   elseif nl==2
     uxx(1)=((-7.)*u(  1)...
            +( 8.)*u(  2)...
            +(-1.)*u(  3))/(2.*dx^2)...
            +(-6.)*ux( 1) /(2.*dx);
      end
%
%  uxx at the right boundary, without ux
   if nu==1
     uxx(n)=(( 2.)*u(n  )...
            +(-5.)*u(n-1)...
            +( 4.)*u(n-2)...
            +(-1.)*u(n-3))/(dx^2);
%
%  uxx at the right boundary, including ux
   elseif nu==2
     uxx(n)=((-7.)*u(n  )...
            +( 8.)*u(n-1)...
            +(-1.)*u(n-2))/(2.*dx^2)...
            +( 6.)*ux(n ) /(2.*dx);
   end
%
%  uxx at the interior grid points
   for i=2:n-1
     uxx(i)=(u(i+1)-2.*u(i)+u(i-1))/dx^2;
   end
%