function dy=der(y,x,n);
 
% dy = der(y,x,n)
%
% DER calculates the n'th difference of y :  dy=y/x 
% using a third order scheme for n=1, fourth order for n=2 and
% fifth order for n=3
%       y       matrix to be derived (data sets are columns)
%	x	x-axis for y, has to be equidistant
%       n       derivative order n=1,2 or 3 (default is 1)
%       dy       derived matrix
 
%                                                       ThDdW  5/89
%changed by Christiane Schroeder 25.4.2003: included x, so that also
%for dx~=1 derivative is correctly calculated		 
 
[m,n2]=size(y);
if m==1, y=y.'; m=n2; n2=1; end
if nargin<2, x=1:m; end
if nargin<3, n=1; end                   % default is first derivative

dx=x(2)-x(1);


if n==1,
        edge = 3*y([1,m],:) - 4*y([2,m-1],:) + y([3,m-2],:);    % y' at edge
        dy = [-edge(1,:);  y(3:m,:)-y(1:m-2,:);  edge(2,:)]/2;
	dy=dy/dx;
elseif n==2,
        edge = 2*y([1,m],:) - 5*y([2,m-1],:) + 4*y([3,m-2],:) ...
                            - y([4,m-3],:);                     % y'' at edge
        dy = [-edge(1,:);  y(3:m,:)-2*y(2:m-1,:)+y(1:m-2,:);  edge(2,:)];
	dy=dy./(dx^2);
elseif n==3,
        edge = (5*y([1:2,m-1:m],:)      - 18*y([2:3,m-2:m-1],:) ...
             + 24*y([3:4,m-3:m-2],:)    - 14*y([4:5,m-4:m-3],:) ...
             +  3*y([5:6,m-5:m-4],:));
        dy = [y(5:m,:)-2*y(4:m-1,:)+2*y(2:m-3,:)-y(1:m-4,:)];
        dy = [-edge(1:2,:); dy; edge(3:4,:)]/2;
	dy=dy./(dx^3);
end