function y = euler_improved(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+1), y(:,n) + h*k1);
    y(:, n+1) = y(:,n) + h/2 * ( k1 + k2);
end

end
