function y = nystroem(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;

for n = 1:N,
    k1 = f( t(n), y(:,n));
    k2 = f( t(n) + 2/3 * h, y(:,n) + 2/3 * h * k1);
    k3 = f( t(n) + 2/3 * h, y(:,n) + 2/3 * h * k2);
    y(:, n+1) = y(:,n) + h/8 * ( 2*k1 + 3*k2 + 3*k3 );
end

end
