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);
i=2:n-1;
dydx(i) = (y(i+1)-y(i-1)) ./ (x(i+1)-x(i-1));
% dydx(1) = (y(2)-y(1)) / (x(2)-x(1));            % simpler
% dydx(n) = (y(n)-y(n-1)) / (x(n)-x(n-1));
dydx(1) = (-y(3)+4*y(2)-3*y(1)) / (x(3)-x(1));  % more accurate
dydx(n) = (3*y(n)-4*y(n-1)+y(n-2)) / (x(n)-x(n-2));

end % function derivative
