% dmFirst5Pt.m
%          
%  1st order 5 point approximation (differentiated degree 4 interpolating polynomial)
%   to the 1st derivative on a uniform grid

%       x    the grid
%      dm    differentiation matrix

% Note: test on a unit spaced grid
%     x = linspace(0,10,11)
%    dm = dmFirst5Pt(x) 

%  the weights come from:
%{
>> format rat
>> C = finiteDifferenceWeights(0,[0 1 2 3 4],1)
C =
       1              0              0              0              0       
     -25/12           4             -3              4/3           -1/4     
>> C = finiteDifferenceWeights(1,[0 1 2 3 4],1)
C =
       0              1              0              0              0       
      -1/4           -5/6            3/2           -1/2            1/12    
>> C = finiteDifferenceWeights(2,[0 1 2 3 4],1)
C =
       0              0              1              0              0       
       1/12          -2/3            0              2/3           -1/12    
>> C = finiteDifferenceWeights(3,[0 1 2 3 4],1)
C =
       0              0              0              1              0       
      -1/12           1/2           -3/2            5/6            1/4     
>> C = finiteDifferenceWeights(4,[0 1 2 3 4],1)
C =
       0              0              0              0              1       
       1/4           -4/3            3             -4             25/12  
%}

 function dm = dmFirst5Pt(x) 

     N = length(x);
     h = abs(x(2)-x(1));
     
     up2 = -diag(ones(1,N),2)/12;
     up = 2*diag(ones(1,N),1)/3;
     low = -2*diag(ones(1,N),-1)/3;
     low2 = diag(ones(1,N),-2)/12;
     
   
     dm = diag(0*ones(1,N)) + up(2:end,2:end) + low(2:end,2:end) + up2(3:end,3:end) + low2(3:end,3:end);
     
     dm(1,1) = -25/12;   dm(1,2) = 4;    dm(1,3) = -3;   dm(1,4) = 4/3;  dm(1,5) = -1/4; 
     dm(2,1) = -1/4;   dm(2,2) = -5/6;    dm(2,3) = 3/2;   dm(2,4) = -1/2;  dm(2,5) = 1/12;
     
     dm(N-1,N-4) = -1/12;   dm(N-1,N-3) = 1/2;     dm(N-1,N-2) = -3/2;     dm(N-1,N-1) = 5/6;    dm(N-1,N) = 1/4; 
       dm(N,N-4) = 1/4;       dm(N,N-3) = -4/3;      dm(N,N-2) = 3;          dm(N,N-1) = -4;       dm(N,N) = 25/12; 
     
     
     dm = dm./h;