function d = diffes(x,y)

% diffes(x,y)
% takes smoother numerical derivative of y with respect to x 
% by fitting quadratics to groups of 3 adjacent points and 
% returns a vector EQUAL IN LENGTH to x and y (unlike diff)

n = length(y);
d = zeros(size(y));
d(1) = (-3*y(1)+4*y(2)-y(3))/(x(3)-x(1));
d(2) = (-3*y(2)+4*y(3)-y(4))/(x(4)-x(2));
for i = 3:n-2,
  d(i) = (-y(i+2)+8*y(i+1)-8*y(i-1)+y(i-2))/(12*(x(2)-x(1)));
end
d(n-1) = (3*y(n-1)-4*y(n-2)+y(n-3))/(x(n-1)-x(n-3));
d(n) = (3*y(n)-4*y(n-1)+y(n-2))/(x(n)-x(n-2));
