function [T , Y] = myeuler(f,tspan,y0,n)
% Solves y' = f(t,y) with initial condition y(a) = y0
% on the interval [a,b] using n steps of Euler's method.
% y must be a column vector and f must produce column vectors
% The dimension of y0 must match that of y and f
% T will be a n+1 column vector contianing the times
% 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
    y = y + h*f(t,y);
    t = a + i*h;
    T = [T; t];
    Y = [Y; y'];
end
plot(T,Y)