function d = diffe(x,y)
% diffe(x,y)
% takes numerical derivative of y with respect to x according
% to dy/dx[i] = (y[i+1]-y[i-1])/(x[i+1]-x[i-1]) and returns a
% vector EQUAL IN LENGTH to x and y (unlike the function diff)
n = length(y);
d = zeros(size(y));
d(1) = (-3*y(1)+4*y(2)-y(3))/(x(3)-x(1));
for i = 2:n-1,
d(i) = (y(i+1)-y(i-1))/(x(i+1)-x(i-1));
end
d(n) = (3*y(n)-4*y(n-1)+y(n-2))/(x(n)-x(n-2));