function [T , Y] = mymodeuler(f,tspan,y0,n)
% Solves y' = f(t,y) with initial condition y(a) = y0
% on the interval [a,b] using n steps of the modified Euler's method.
% y and f may be column vectors
% The dimension of y0 must match that of y and f
% T will be a n+1 column vector
% Y will be a (n+1) by d matrix where d is the length of y
a = tspan(1);
b = tspan(2);
h = (b-a)/n; % step size
t = a; y = y0; % t and y are the current variables
T = a; Y = y0'; % T and Y will record all steps
for i = 1:n
k1 = h*f(t,y);
k2 = h*f(t+h,y+k1);
y = y + .5*(k1+k2);
t = a + i*h;
T = [T; t];
Y = [Y; y'];
end
plot(T,Y)
