% Copyright 2008 by Jonathan A. Cox, joncox@mit.edu. All rights reserved.
% Unauthorized reproduction or distribution is prohibited.

% 1st Order Derivate: Central difference, O(h^2)

function d = cdiff1(y,h)

%d = zeros(size(y));

%for k=2:length(y)-1
%    d(k) = (y(k+1)-y(k-1))/(2*h);
%end

% This is only slightly faster than the for loop
%d = (circshift(y,[1 -1]) - circshift(y,[1 1]))/(2*h);

% This is by far the fastest method of taking the derivative
d = (y([2:end 1]) - y([end 1:end-1]))/(2*h);

d(1) = d(2) + (d(2)-d(3));
d(end) = d(end-1) + (d(end-1)-d(end-2));