function [t,y] = eulertd2(f,t0,tf,y0,n)
% eulertd2.m - forward euler for systems of ODEs
% T. DeLillo, Math 451, Fall 2001
  global wsq
  wsq = 1.56; % specify value of w^2 = g/l;
  h = (tf-t0)/n;
  t = (t0:h:tf);
  y = y0;
% forward euler loop
  for i = 1:n
     y(:,i+1)=y(:,i)+h*feval(f,t(i),y(:,i));
  end