function dydx=derivative(x,y)

% Given a function y as a vector of values at points x, 
% it returns its derivative at the same points. 
% Uses a simple finite-difference approximation:
%   dydx(i) = (y(i+1)-y(i-1)) / (x(i+1)-x(i-1))   for 1<i<n
%   dydx(1) = (y(2)-y(1)) / (x(2)-x(1))
%   dydx(n) = (y(n)-y(n-1)) / (x(n)-x(n-1))
% Input:
%   x    : vector of x points
%   y    : vector of y values at given x points
% Output:
%   dydx : vector of dy/dx at given x points 
% Usage example:
%    x = 0:pi/100:2*pi
%    y = sin(x)
%    z = derivative(x,y)
% J.M.Soler. Sep.2010

n=numel(x);
% for i=1:n
%     if i==1
%         dydx(i) = (y(2)-y(1)) / (x(2)-x(1));            % simpler
% %        dydx(i) = (-y(3)+4*y(2)-3*y(1)) / (x(3)-x(1));  % more accurate
%     elseif i==n
%         dydx(i) = (y(n)-y(n-1)) / (x(n)-x(n-1));
% %        dydx(i) = (3*y(n)-4*y(n-1)+y(n-2)) / (x(n)-x(n-2));
%     else
%         dydx(i) = (y(i+1)-y(i-1)) / (x(i+1)-x(i-1));
%     end
% end

% A faster alternative (and more accurate for non-uniform mesh)
i = 2:n-1;                             % all points except first and last
dxdi(i) = (x(i+1)-x(i-1))/2;           % dx(i)/di
dydi(i) = (y(i+1)-y(i-1))/2;           % dy(i)/di
dxdi(1) = (-x(3)+4*x(2)-3*x(1))/2;     % assumes quadratic x(i)
dydi(1) = (-y(3)+4*y(2)-3*y(1))/2;
dxdi(n) = (3*x(n)-4*x(n-1)+x(n-2))/2;
dydi(n) = (3*y(n)-4*y(n-1)+y(n-2))/2;
dydx = dydi./dxdi;                     % dy(i)/dx = (dy(i)/di)/(dx(i)/di)

end % function derivative
