  function [t,y] = eulertd3(f,t0,tf,y0,n)
% eulertd3.m - forward euler for systems of ODEs
% for pendulum with animated graphics in (x,y)
% T. DeLillo, Math 451, Fall 2001
  clf
  global wsq b
  l=1;
  b=0;
  wsq = 1.56/l; % specify value of w^2 = g/l;
  h = (tf-t0)/n;
  t = (t0:h:tf);
  y = y0;   % initial conditions
% forward euler loop
     axis([-l-.5 l+.5 -l-.5 l+.5]);
     axis('square');
     hbead=line(t(1),y(1,1),'marker','o','markersize',8,'erase','xor');
     set(hbead,'xdata',l*sin(y(1,1)),'ydata',-l*cos(y(1,1)));
     drawnow
     pause(.1)
  for i = 1:n
     y(:,i+1)=y(:,i)+h*feval(f,t(i),y(:,i));
     set(hbead,'xdata',l*sin(y(1,i+1)),'ydata',-l*cos(y(1,i+1)));
     drawnow
     pause(.1)
  end 