function [dudt] = accel(u,dt);
% ACCEL - Calculate acceleration
% [dudt] = accel(u,dt);
%
% First derivative using 5 pt centered difference
% Ignore two points on each end of returned array.

% Chris Sherwood, USGS
n = length(u);
m2 = [u(1);u(1);u(1:n-2)];
m1 = [u(1);u(1:n-1)];
p1 = [u(2:n);u(n)];
p2 = [u(3:n);u(n);u(n)];
% dudt = (p1-m1)./2*dt; % 3 pt centered diff
dudt = (-p2 + 8*p1 - 8*m1 + m2)./12*dt; % 5 pt centered diff