function [xpts,y]=euler(fun,a,b,N,y0)

% EULER  Euler's method 
%        [xpts,y] = euler(fun,a,b,N,y0) returns the numerical solution to
%        y'=f(x,y) for x in [a,b] with initial condition x(a)=y0.
%        fun is the name of a function must accept x and y as inputs
%        N is the number of intervals to use on [a,b]

h=(b-a)/N;
y=zeros(N+1,1); y(1)=y0;
% grid of x-pts x_{n}
xpts=a:h:b;
for i=2:N+1
    y(i)=y(i-1)+h*feval(fun, xpts(i-1),y(i-1));
end

