function y = euler_modified(h, f, t0, y0, T)
%EULER_EXPLICIT Summary of this function goes here
%   Detailed explanation goes here

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

for n = 1:N,
    t_mid = t(n) + 0.5 * h;
    y_mid = y(:,n) + 0.5 * h * f( t(n), y(:,n)); 
    y(:, n+1) = y(:,n) + h * f( t_mid, y_mid);
end

end

