function yi = myDevalInt2(sol, ix, 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

if ix == 1
    yi = zeros(length(inds),1);
    return
end

% The matrix to transform coefficients to powers of s

A = [  1   1/2   1/3   1/4    1/5
       0   1/2   1/2 11/24   5/12
       0     0   1/6   1/4   7/24
       0     0     0  1/24   1/12
       0     0     0     0  1/120]';

k    = sol.idata.kvec(ix+1);
K    = (1:k)';

D    = sol.idata.dif3d(:,:,ix+1);
C    = D(inds,K)*A(K,K);

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

s    = (t - tnew)/h;

scale = [1/2;1/3;1/4;1/5;1/6];
yi = ynew*h + C*(scale(1:k).*s(ones(k,1)).*cumprod(s(ones(k,1))))*h;