  function [t,y] = eulertd1(f,t0,tf,y0,n)
% eulerd.m - forward euler
  h = (tf-t0)/n;
  t = (t0:h:tf)';
  y(1)=y0;
% forward euler loop
  for i = 1:n
     y(i+1)=y(i)+h*feval(f,t(i),y(i));
  end 