% File: dss050.m
%
  function [uxx]=dss050(xl,xu,n,u,ux,nl,nu)
%
%  Function dss050 computes a tenth-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 tenth-order, twelve-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);
%
%  1/dx**2 for subsequent use
   rdxs=1.0/dx^2;
%
%  uxx at the left boundary
%
%     Without ux 
         if nl==1
         uxx(1)=rdxs*...
               (  7.561626984126982*u(1)...
                -44.437301587301590*u(2)...
               +138.593253968253980*u(3)...
               -295.519841269841320*u(4)...
               +457.029761904761930*u(5)...
               -521.113333333333460*u(6)...
               +439.394444444444330*u(7)...
               -271.261904761904700*u(8)...
               +119.413690476190470*u(9)...
                -35.551587301587304*u(10)...
                 +6.423730158730158*u(11)...
                 -0.532539682539683*u(12));
%
%     With ux
         elseif nl==2
         uxx(1)=rdxs*...
               (-10.128622763920379*u(1)...
                +20.000000000000000*u(2)...
                -22.499999999999972*u(3)...
                +26.666666666666572*u(4)...
                -26.249999999999943*u(5)...
                +20.159999999999968*u(6)...
                -11.666666666666686*u(7)...
                 +4.897959183673493*u(8)...
                 -1.406250000000000*u(9)... 
                 +0.246913580246911*u(10)... 
                 -0.020000000000001*u(11)... 
                 -5.857936507936508*ux(1)*dx);
         end
%
%  uxx at the right boundary
%
%     Without ux 
         if nu==1
         uxx(n)=rdxs*...
               ( -0.532539682539683*u(n-11)...
                 +6.423730158730161*u(n-10)...
                -35.551587301587304*u(n-9)...
               +119.413690476190480*u(n-8)...
               -271.261904761904760*u(n-7)...
               +439.394444444444390*u(n-6)...
               -521.113333333333340*u(n-5)...
               +457.029761904761930*u(n-4)...
               -295.519841269841270*u(n-3)...
               +138.593253968253980*u(n-2)...
                -44.437301587301590*u(n-1)...
                 +7.561626984126984*u(n  ));
%
%     With ux 
         elseif nu==2
         uxx(n)=rdxs*...
               ( -0.019999999999999*u(n-10)...
                 +0.246913580246904*u(n-9)...
                 -1.406249999999986*u(n-8)...
                 +4.897959183673493*u(n-7)...
                -11.666666666666686*u(n-6)...
                +20.159999999999968*u(n-5)...
                -26.249999999999943*u(n-4)...
                +26.666666666666629*u(n-3)...
                -22.499999999999972*u(n-2)...
                +20.000000000000000*u(n-1)...
                -10.128622763920383*u(n  )...
                 +5.857936507936508*ux(n  )*dx);
         end
%
%  uxx at the interior grid points
%
%     i = 2 
         uxx(2)=rdxs*...
               (  0.532539682539683*u(1)...
                 +1.171150793650794*u(2)...
                 -9.289682539682540*u(3)...
                +21.434523809523807*u(4)...
                -31.912698412698408*u(5)...
                +35.258333333333333*u(6)...
                -29.046666666666663*u(7)...
                +17.623015873015873*u(8)...
                 -7.654761904761902*u(9)...
                 +2.254960317460317*u(10)...
                 -0.403968253968254*u(11)...
                 +0.033253968253968*u(12));
%
%     i = 3
         uxx(3)=rdxs*...
               ( -0.033253968253968*u(1)...
                 +0.931587301587302*u(2)...
                 -1.023611111111112*u(3)...
                 -1.973809523809523*u(4)...
                 +4.973809523809524*u(5)...
                 -5.575555555555556*u(6)...
                 +4.531666666666666*u(7)...
                 -2.709523809523809*u(8)...
                 +1.162301587301588*u(9)...
                 -0.338888888888889*u(10)...
                 +0.060198412698413*u(11)...
                 -0.004920634920635*u(12));
%
%     i = 4
         uxx(4)=rdxs*...
               (  0.004920634920635*u(1)...
                 -0.092301587301587*u(2)...
                 +1.256349206349206*u(3)...
                 -2.106150793650793*u(4)...
                 +0.461904761904762*u(5)...
                 +1.076666666666667*u(6)...
                 -1.028888888888889*u(7)...
                 +0.634523809523810*u(8)...
                 -0.273809523809524*u(9)...
                 +0.079761904761905*u(10)...
                 -0.014126984126984*u(11)...
                 +0.001150793650794*u(12));
%
%     i = 5
         uxx(5)=rdxs*...
               ( -0.001150793650794*u(1)...
                 +0.018730158730159*u(2)...
                 -0.168253968253968*u(3)...
                 +1.509523809523809*u(4)...
                 -2.675793650793650*u(5)...
                 +1.373333333333333*u(6)...
                 +0.013333333333334*u(7)...
                 -0.117460317460318*u(8)...
                 +0.064880952380952*u(9)...
                 -0.020634920634921*u(10)...
                 +0.003809523809524*u(11)...
                 -0.000317460317460*u(12));
%
%     i = n-4
         uxx(n-4)=rdxs*...
                (-0.000317460317460*u(n-11)...
                 +0.003809523809524*u(n-10)...
                 -0.020634920634921*u(n-9)...
                 +0.064880952380952*u(n-8)...
                 -0.117460317460317*u(n-7)...
                 +0.013333333333333*u(n-6)...
                 +1.373333333333333*u(n-5)...
                 -2.675793650793652*u(n-4)...
                 +1.509523809523810*u(n-3)...
                 -0.168253968253968*u(n-2)...
                 +0.018730158730159*u(n-1)...
                 -0.001150793650794*u(n  ));
%
%     i = n-3
         uxx(n-3)=rdxs*...
                ( 0.001150793650794*u(n-11)...
                 -0.014126984126984*u(n-10)...
                 +0.079761904761905*u(n-9)...
                 -0.273809523809524*u(n-8)...
                 +0.634523809523809*u(n-7)...
                 -1.028888888888889*u(n-6)...
                 +1.076666666666667*u(n-5)...
                 +0.461904761904762*u(n-4)...
                 -2.106150793650793*u(n-3)...
                 +1.256349206349206*u(n-2)...
                 -0.092301587301587*u(n-1)...
                 +0.004920634920635*u(n  ));
%
%     i = n-2
         uxx(n-2)=rdxs*...
                (-0.004920634920635*u(n-11)...
                 +0.060198412698413*u(n-10)...
                 -0.338888888888889*u(n-9)...
                 +1.162301587301588*u(n-8)...
                 -2.709523809523809*u(n-7)...
                 +4.531666666666666*u(n-6)...
                 -5.575555555555556*u(n-5)...
                 +4.973809523809526*u(n-4)...
                 -1.973809523809526*u(n-3)...
                 -1.023611111111110*u(n-2)...
                 +0.931587301587302*u(n-1)...
                 -0.033253968253968*u(n  ));
%
%     i = n-1
         uxx(n-1)=rdxs*...
                ( 0.033253968253968*u(n-11)...
                 -0.403968253968254*u(n-10)...
                 +2.254960317460318*u(n-9)...
                 -7.654761904761904*u(n-8)...
                +17.623015873015873*u(n-7)...
                -29.046666666666670*u(n-6)...
                +35.258333333333340*u(n-5)...
                -31.912698412698411*u(n-4)...
                +21.434523809523810*u(n-3)...
                 -9.289682539682541*u(n-2)...
                 +1.171150793650794*u(n-1)...
                 +0.532539682539683*u(n  ));
%
%     i = 6, 7,..., n-5
         for i=6:n-5
         uxx(i)=rdxs*...
                ( 0.000317460317460*u(i-5)...
                 -0.004960317460317*u(i-4)...
                 +0.039682539682540*u(i-3)...
                 -0.238095238095238*u(i-2)...
                 +1.666666666666667*u(i-1)...
                 -2.927222222222222*u(i  )...
                 +1.666666666666667*u(i+1)...
                 -0.238095238095238*u(i+2)...
                 +0.039682539682540*u(i+3)...
                 -0.004960317460317*u(i+4)...
                 +0.000317460317460*u(i+5));
         end

