function [y,t] = modified_euler(f,t0,y0,T,n)
% by Shelvean Kapita
% modified Euler's method to solve ODE
% y' = f(t,y)
% t0 = initial time
% y0 = initial value
% T  = final time
% n = number of time steps
h = (T-t0)/n; % time step
i = 1;
t = linspace(t0,T,n+1);
y = zeros(size(t));
y(1)=y0;
while (i<=n)
    f1 = f(t(i),y(i));
    y1 = y(i)+h*f(t(i),y(i));
    f2 = f(t(i+1),y1);
    y(i+1) = y(i)+0.5*h*(f1+f2);
    i = i+1;
end