function [d,z,finitedifference] = residual( t, y, fun, refin )
%
% Compute z'(t) - fun(t,z(t))
% where t, y is the ode45 solution to y' = fun(t,y)
% with the "refine" option refin (default 4)
%
if nargin==3,
    ref = 4;
else
    ref = refin;
end;
top = [ -25/12, 4, -3, 4/3, -1/4,
        -1/4, -5/6, 3/2, -1/2, 1/12];
middle = [1/12, -2/3, 0, 2/3, -1/12];
bottom = [-1/12, 1/2, -3/2, 5/6, 1/4,
           1/4, -4/3, 3, -4, 25/12];
finitedifference = zeros(ref+1,ref+1);
finitedifference(1:2,1:5) = top;
for i=3:ref-1,
    finitedifference(i,i-2:i+2) = middle;
end;
finitedifference(ref:ref+1,ref-3:ref+1) = bottom;
[m,n] = size(y); % m == ell
z = zeros(m,n);
k = (m-1)/ref;
for i=1:k,
    j = ref*(i-1);
    yblock = y(1+j:ref+1+j,:);
    h = t(2+j)-t(1+j);
    zblock = finitedifference*yblock/h;
    z(1+j:ref+1+j,:) = zblock; % wipes out last row
end;
d = zeros(m,n);
for i=1:m,
    d(i,:) = z(i,:) - feval(fun,t(i),y(i,:));
end;