function [ylist,tlist] = euler(y0, t0, m, vf, h, N)
% solves dy/dt = f(y,t) in R^m starting from y(t0)=y0, using N steps of size h
% the initial condition should be an m-dimensional column vector
% vf(y,t) should evaluate the vector field f at the point (y,t)
% results are returned in tlist and ylist, the latter in the form of a matrix, each
% of whose columns represent the numerical solution at a timestep

tlist = zeros(1,N+1);
ylist = zeros(m, N+1);
tlist(1) = t0;
ylist(:,1) = y0;

for nn=1:N,
	tlist(nn+1) = t0 + nn*h;
	ylist(:,nn+1) = ylist(:,nn) + h*feval(vf, tlist(:,nn), ylist(:,nn));
end
