function [ t, resid, yhat, yphat ] = ODEResidualMass( sol, M, f, nRefine )
%Compute the residual of the solution sol in the de My' = f(t,y)
%   ODEResidualMass( sol, M, f, nRefine )
    if nargin < 4,
        nRefine = 4;
    end
    t = RefineMesh( sol.x, nRefine );
    nump = length(t);
    
    [yhat,yphat] = deval( sol, t );
    resid = zeros( size(yphat) );
    yphat = M*yphat;
    for i=1:nump,
        resid(:,i) = yphat(:,i) - f(t(i), yhat(:,i));
    end
    
end

