function [ X,Y ] = MyEuler( odefun, span, y0, n )
%Implements Euler ODE solver for
% dy/dx = odefun(x,y) with y(a) = y0 (a column vec)

% Interval
a = span(1);
b = span(2);

x = a;
X = a;
y = y0;
Y = y0';   % vector equation

h = (b-a)/n;

for i = 1:n
    y = y + h * odefun(x,y); 
    % y2 = y + h * odefun(x,y);
    % y = y + h * 1/2 * (odefun(x,y) + odefun(x+h,y2))
    x = a + i*h;
    X = [X; x];
    Y = [Y; y'];
end

% g = inline('[y(2); -y(1)]', 'x','y')




end