% File: dss048.m
%
  function [uxx]=dss048(xl,xu,n,u,ux,nl,nu)
%
%  Function dss048 computes a eighth-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 eighth-order, ten-point approximations for a
%  second derivative 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.
%
%  Weighting coefficients for finite difference approximations, computed
%  by an algorithm of B. Fornberg (1,2), are used in the following 
%  approximations.  The combination of these coefficients to give the 
%  second derivative approximations at the boundaries in x, including 
%  the first derivative, was suggested by Professor Fornberg.
%
%  (1)  Fornberg, B., Fast Generation of Weights in Finite Difference
%       Formulas, In Recent Developments In Numerical Methods And
%       Software for ODEs/DAEs/PDEs, G. Byrne and W. E. Schiesser (eds), 
%       World Scientific, River Edge, NJ, 1992
%
%  (2)  Fornberg, B., Calculation of Weights In Finite Difference
%       Formulas, SIAM Review, Vol. 40, No. 3, Pp 685-691, September, 
%       1999
%
%  Grid spacing
   dx=(xu-xl)/(n-1);
   rdxs=1.0/dx^2;
%
%  uxx at the left boundary
%
%     Without ux 
         if nl==1
         uxx(1)=rdxs*...
                ( 6.463293650793652*u(1)...
                -32.921428571428578*u(2)...
                +83.842857142857156*u(3)...
               -139.755555555555555*u(4)...
               +162.474999999999999*u(5)...
               -132.500000000000000*u(6)...
                +74.544444444444444*u(7)...
                -27.628571428571430*u(8)...
                 +6.083928571428572*u(9)...
                 -0.603968253968254*u(10));
%
%     With ux
         elseif nl==2
         uxx(1)=rdxs*...
               ( -8.914169501133784*u(1)...
                +16.000000000000000*u(2)...
                -14.000000000000000*u(3)...
                +12.444444444444444*u(4)...
                 -8.750000000000000*u(5)...
                 +4.480000000000000*u(6)...
                 -1.555555555555555*u(7)...
                 +0.326530612244898*u(8)...
                 -0.031250000000000*u(9)... 
                 -5.435714285714286*ux(1)*dx);
         end
%
%  uxx at the right boundary
%
%     Without ux 
         if nu==1
         uxx(n)=rdxs*...
               ( -0.603968253968254*u(n-9)...
                 +6.083928571428572*u(n-8)...
                -27.628571428571430*u(n-7)...
                +74.544444444444444*u(n-6)...
               -132.500000000000000*u(n-5)...
               +162.474999999999999*u(n-4)...
               -139.755555555555555*u(n-3)...
                +83.842857142857156*u(n-2)...
                -32.921428571428578*u(n-1)...
                 +6.463293650793652*u(n  ));
%
%     With ux 
         elseif nu==2
         uxx(n)=rdxs*...
                (-0.031250000000000*u(n-8)...
                 +0.326530612244905*u(n-7)...
                 -1.555555555555555*u(n-6)...
                 +4.480000000000000*u(n-5)...
                 -8.750000000000000*u(n-4)...
                +12.444444444444444*u(n-3)...
                -14.000000000000000*u(n-2)...
                +16.000000000000000*u(n-1)...
                 -8.914169501133790*u(n  )...
                 +5.435714285714287*ux(n  )*dx);
         end
%
%  uxx at the interior grid points
%
%     i = 2 
         uxx(2)=rdxs*...
                 (0.603968253968254*u(1)...
                 +0.423611111111111*u(2)...
                 -5.742857142857143*u(3)...
                +11.366666666666665*u(4)...
                -12.922222222222221*u(5)...
                +10.274999999999999*u(6)...
                 -5.666666666666665*u(7)...
                 +2.068253968253968*u(8)...
                 -0.450000000000000*u(9)...
                 +0.044246031746032*u(10));
%
%     i = 3
         uxx(3)=rdxs*...
                (-0.044246031746032*u(1)...
                 +1.046428571428571*u(2)...
                 -1.567460317460318*u(3)...
                 -0.433333333333333*u(4)...
                 +2.075000000000000*u(5)...
                 -1.772222222222222*u(6)...
                 +0.983333333333333*u(7)...
                 -0.357142857142857*u(8)...
                 +0.077182539682540*u(9)...
                 -0.007539682539683*u(10));
%
%     i = 4
         uxx(4)=rdxs*...
                 (0.007539682539683*u(1)...
                 -0.119642857142857*u(2)...
                 +1.385714285714286*u(3)...
                 -2.472222222222222*u(4)...
                 +1.150000000000000*u(5)...
                 +0.175000000000000*u(6)...
                 -0.188888888888889*u(7)...
                 +0.078571428571429*u(8)...
                 -0.017857142857143*u(9)...
                 +0.001785714285714*u(10));
%
%     i = n-3
         uxx(n-3)=rdxs*...
                 (0.001785714285714*u(n-9)...
                 -0.017857142857143*u(n-8)...
                 +0.078571428571429*u(n-7)...
                 -0.188888888888889*u(n-6)...
                 +0.175000000000000*u(n-5)...
                 +1.150000000000000*u(n-4)...
                 -2.472222222222224*u(n-3)...
                 +1.385714285714286*u(n-2)...
                 -0.119642857142857*u(n-1)...
                 +0.007539682539683*u(n  ));
%
%     i = n-2 
         uxx(n-2)=rdxs*...
                (-0.007539682539683*u(n-9)...
                 +0.077182539682540*u(n-8)...
                 -0.357142857142857*u(n-7)...
                 +0.983333333333333*u(n-6)...
                 -1.772222222222222*u(n-5)...
                 +2.075000000000001*u(n-4)...
                 -0.433333333333335*u(n-3)...
                 -1.567460317460318*u(n-2)...
                 +1.046428571428572*u(n-1)...
                 -0.044246031746032*u(n  ));
%
%     i = n-1 
         uxx(n-1)=rdxs*...
                 (0.044246031746032*u(n-9)...
                 -0.450000000000000*u(n-8)...
                 +2.068253968253968*u(n-7)...
                 -5.666666666666667*u(n-6)...
                +10.275000000000000*u(n-5)...
                -12.922222222222222*u(n-4)...
                +11.366666666666667*u(n-3)...
                 -5.742857142857142*u(n-2)...
                 +0.423611111111112*u(n-1)...
                 +0.603968253968254*u(n  ));
%
%     i = 5, 6,..., n-4
         for i=5:n-4
         uxx(i)=rdxs*...
                (-0.001785714285714*u(i-4)...
                 +0.025396825396825*u(i-3)...
                 -0.200000000000000*u(i-2)...
                 +1.600000000000000*u(i-1)...
                 -2.847222222222221*u(i  )...
                 +1.600000000000000*u(i+1)...
                 -0.200000000000000*u(i+2)...
                 +0.025396825396825*u(i+3)...
                 -0.001785714285714*u(i+4));
         end
