function y = TwoStep(h, f, t0, y0, T)

t = t0:h:T;
N = length(t)-1;
d = size(y0,1);
y = zeros(d, N+1);
y(:,1) = y0;
options=[];

%options=odeset('RelTol',1e-14,'AbsTol',1e-14);
[~,yh] = ode45(f,t(1:2),y0,options);
y(:,2) = yh(end,:)';

for n = 1:N-1,
    y(:, n+2) = y(:,n) + 2 * h * f( t(n+1), y(:,n+1));
end

end
