deriv

 DERIV Column-wise derivative estimation.
    Computes 5-point discrete derivative estimates for each column
    of the input matrix U.
    V = DERIV(U)   Returns a matrix V of the same size as U each
    element containing the estimates of the column-wise derivative of
    a function sampled at unit intervals in U.

    DERIV Uses 5 points both in the interior and at the edges.
    If size of the matrix U is less than 5 use simpler 3-point
    or 2-point estimate.

    See also DIFF, GRADIENT, DEL2

0001 function  V = deriv(U)
0002 
0003 % DERIV Column-wise derivative estimation.
0004 %    Computes 5-point discrete derivative estimates for each column
0005 %    of the input matrix U.
0006 %    V = DERIV(U)   Returns a matrix V of the same size as U each
0007 %    element containing the estimates of the column-wise derivative of
0008 %    a function sampled at unit intervals in U.
0009 %
0010 %    DERIV Uses 5 points both in the interior and at the edges.
0011 %    If size of the matrix U is less than 5 use simpler 3-point
0012 %    or 2-point estimate.
0013 %
0014 %    See also DIFF, GRADIENT, DEL2
0015 
0016 %  Kirill K. Pankratov,  kirill@plume.mit.edu
0017 %  March 22, 1994
0018 
0019  % 5-point coefficients for derivatives ..................................
0020  % Edges (1st, 2nd pts):
0021 Cde = [-25/12 4 -3 4/3 -1/4; -1/4 -5/6 3/2 -1/2 1/12];
0022  % Interior
0023 Cd = [1/12 -2/3 0 2/3 -1/12];
0024 
0025  % Determine the size of the input and make column if vector .............
0026 ist = 0;
0027 ly = size(U,1);
0028 if ly==1, ist = 1; U = U(:); ly = length(U); end
0029 
0030  % If only 2 points - simple difference ..................................
0031 if ly==2
0032   V(1,:) = U(2,:)-U(1,:);
0033   V(2,:) = V(1,:);
0034   if ist, V = V'; end     % Transpose output if necessary
0035   return
0036 end
0037 
0038  % Now if more than 2 points - more complicated procedure ...............
0039 if ly<5         % If less than 5 points
0040 
0041   V(2:ly-1,:) = (U(3:ly,:)-U(1:ly-2,:))/2;        % First
0042   V(1,:) = 2*U(2,:)-1.5*U(1,:)-.5*U(3,:);         % Last
0043   V(ly,:) = 1.5*U(ly,:)-2*U(ly-1,:)+.5*U(ly-2,:); % Interior
0044 
0045 else            % If 5 points or more - regular, more accurate estimation
0046 
0047    % Edges ............
0048   V(1:2,:) = Cde*U(1:5,:);
0049   V([ly ly-1],:) = -fliplr(Cde)*U(ly-4:ly,:);
0050 
0051    % Interior .........
0052   V(3:ly-2,:) = Cd(1)*U(1:ly-4,:)+Cd(2)*U(2:ly-3,:)+Cd(3)*U(3:ly-2,:);
0053   V(3:ly-2,:) = V(3:ly-2,:)+Cd(4)*U(4:ly-1,:)+Cd(5)*U(5:ly,:);
0054 
0055 end
0056 
0057 if ist, V = V'; end     % Transpose output if necessary