function yi = myDeval(sol, ti, inds)
% This code returns the same answer as ode15s, for a single scalar
% timepoint.

if nargin < 3 || isempty(inds)
    inds = 1:size(sol.y, 1);
end

ix   = find(ti-sol.x >= 0, 1, 'last');

if ti == sol.x(ix)
    yi = sol.y(inds,ix);
    return
end

dif  = sol.idata.dif3d(:,:,ix+1);
k = sol.idata.kvec(ix+1);

K = (1:k)';

t    = sol.x(ix);
tnew = sol.x(ix+1);
ynew = sol.y(inds,ix+1);
h    = tnew-t;

dt   = (ti - tnew)/h;

kI = K;
yi = ynew + dif(inds,1:k)*cumprod((dt(ones(k,1))+kI-1)./kI);