%  File: dss010.m
%
   function [ux]=dss010(xl,xu,n,u)
%
%  Function dss010 computes the first derivative, u , of a
%                                                 x
%  variable u over the spatial domain xl le x le xu from classical
%  eleven-point, tenth-order finite difference approximations
%
%  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)
%
%  The mathematical details of the finite difference approxima-
%  tions are given in routines dss002, dss004, dss006 and
%  dss008, and can be summarized by the following Bickley matrix
%  for n = 10, m = 1, p = 0 to 10 (Bickley, W. G., Formulae for
%  Numerical Differentiation, Math. Gaz., vol. 25, 1941)
%
%      -10626840            36288000           -81648000
%      145152000          -190512000           182891520
%
%        -362880            -6636960            16329600
%      -21772800            25401600           -22861440
%
%          40320             -806400            -4419360
%        9676800            -8467200             6773760
%
%         -10080              151200            -1360800
%       -2756160             6350400            -3810240
%
%           4320              -57600              388800
%       -2073600            -1330560             4354560
%
%          -2880               36000             -216000
%         864000            -3024000                   0
%
%           2880              -34560              194400
%        -691200             1814400            -4354560
%
%          -4320               50400             -272160
%         907200            -2116800             3810240
%
%          10080             -115200              604800
%       -1935360             4233600            -6773760
%
%         -40320              453600            -2332800
%        7257600           -15240960            22861440
%
%         362880            -4032000            20412000
%      -62208000           127008000          -182891520
%
%  Each entry in this table should be multiplied by 1/10 to
%  obtain the final weighting coefficients.  From this Bickley
%  matrix, the finite difference approximation of the first
%  derivative can be programmed for each of the grid points 1, 2,
%  3, 4, 5,..., i,..., n-4, n-3, n-2, n-1 and n (taking into
%  account the symmetry properties of the matrix).
%
%  Compute the spatial increment
   dx=(xu-xl)/(n-1);
   r10fdx=1./(3628800.*dx);
   nm5=n-5;
%
%  Grid point 1
   ux(  1)=r10fdx*...
     (-10628640.      *u(   1)...
      +36288000.      *u(   2)...
      -81648000.      *u(   3)...
      +145152000.     *u(   4)...
      -190512000.     *u(   5)...
      +182891520.     *u(   6)...
      -127008000.     *u(   7)...
      +62208000.      *u(   8)...
      -20412000.      *u(   9)...
      +4032000.       *u(  10)...
      -362880.        *u(  11));
%
%  Grid point 2
   ux(  2)=r10fdx*...
     (-362880.        *u(   1)...
      -6636960.       *u(   2)...
      +16329600.      *u(   3)...
      -21772800.      *u(   4)...
      +25401600.      *u(   5)...
      -22861440.      *u(   6)...
      +15240960.      *u(   7)...
      -7257600.       *u(   8)...
      +2332800.       *u(   9)...
      -453600.        *u(  10)...
      +40320.         *u(  11));
%
%  Grid point 3
   ux(  3)=r10fdx*...
     (+40320.         *u(   1)...
      -806400.        *u(   2)...
      -4419360.       *u(   3)...
      +9676800.       *u(   4)...
      -8467200.       *u(   5)...
      +6773760.       *u(   6)...
      -4233600.       *u(   7)...
      +1935360.       *u(   8)...
      -604800.        *u(   9)...
      +115200.        *u(  10)...
      -10080.         *u(  11));
%
%  Grid point 4
   ux(  4)=r10fdx*...
     (-10080.         *u(   1)...
      +151200.        *u(   2)...
      -1360800.       *u(   3)...
      -2756160.       *u(   4)...
      +6350400.       *u(   5)...
      -3810240.       *u(   6)...
      +2116800.       *u(   7)...
      -907200.        *u(   8)...
      +272160.        *u(   9)...
      -50400.         *u(  10)...
      +4320.          *u(  11));
%
%  Grid point 5
   ux(  5)=r10fdx*...
     (+4320.          *u(   1)...
      -57600.         *u(   2)...
      +388800.        *u(   3)...
      -2073600.       *u(   4)...
      -1330560.       *u(   5)...
      +4354560.       *u(   6)...
      -1814400.       *u(   7)...
      +691200.        *u(   8)...
      -194400.        *u(   9)...
      +34560.         *u(  10)...
      -2880.          *u(  11));
%
%  Grid point i, i ne 1 to 5, n-4 to n
   for i=6:nm5
     ux(  i)=r10fdx*...
     (-2880.          *u( i-5)...
      +36000.         *u( i-4)...
      -216000.        *u( i-3)...
      +864000.        *u( i-2)...
      -3024000.       *u( i-1)...
      +0.             *u( i  )...
      +3024000.       *u( i+1)...
      -864000.        *u( i+2)...
      +216000.        *u( i+3)...
      -36000.         *u( i+4)...
      +2880.          *u( i+5));
   end
%
%  Grid point n-4
   ux(n-4)=r10fdx*...
     (+2880.          *u(n-10)...
      -34560.         *u( n-9)...
      +194400.        *u( n-8)...
      -691200.        *u( n-7)...
      +1814400.       *u( n-6)...
      -4354560.       *u( n-5)...
      +1330560.       *u( n-4)...
      +2073600.       *u( n-3)...
      -388800.        *u( n-2)...
      +57600.         *u( n-1)...
      -4320.          *u( n  ));
%
%  Grid point n-3
   ux(n-3)=r10fdx*...
     (-4320.          *u(n-10)...
      +50400.         *u( n-9)...
      -272160.        *u( n-8)...
      +907200.        *u( n-7)...
      -2116800.       *u( n-6)...
      +3810240.       *u( n-5)...
      -6350400.       *u( n-4)...
      +2756160.       *u( n-3)...
      +1360800.       *u( n-2)...
      -151200.        *u( n-1)...
      +10080.         *u( n  ));
%
%  Grid point n-2
   ux(n-2)=r10fdx*...
     (+10080.         *u(n-10)...
      -115200.        *u( n-9)...
      +604800.        *u( n-8)...
      -1935360.       *u( n-7)...
      +4233600.       *u( n-6)...
      -6773760.       *u( n-5)...
      +8467200.       *u( n-4)...
      -9676800.       *u( n-3)...
      +4419360.       *u( n-2)...
      +806400.        *u( n-1)...
      -40320.         *u( n  ));
%
%  Grid point n-1
   ux(n-1)=r10fdx*...
     (-40320.         *u(n-10)...
      +453600.        *u( n-9)...
      -2332800.       *u( n-8)...
      +7257600.       *u( n-7)...
      -15240960.      *u( n-6)...
      +22861440.      *u( n-5)...
      -25401600.      *u( n-4)...
      +21772800.      *u( n-3)...
      -16329600.      *u( n-2)...
      +6636960.       *u( n-1)...
      +362880.        *u( n  ));
%
%  Grid point n
   ux(n  )=r10fdx*...
     (+362880.        *u(n-10)...
      -4032000.       *u( n-9)...
      +20412000.      *u( n-8)...
      -62208000.      *u( n-7)...
      +127008000.     *u( n-6)...
      -182891520.     *u( n-5)...
      +190512000.     *u( n-4)...
      -145152000.     *u( n-3)...
      +81648000.      *u( n-2)...
      -36288000.      *u( n-1)...
      +10628640.      *u( n  ));
