% Script to numerically solve the pendulum ODE IVP
%
%    u" + sin(u) = 0, u(0) = eps, u'(0) = 0,  0 < t < tmax
%
% for specified a and tmax, plot the solution y(t), and compare
% with the perturbation series from class (truncated after O(eps^3) terms).
% for  eps = 0.5 and 1;

  tmax = 4*pi;
  tspan = [0 tmax];

  epsv = [0.5 1];
  for j = 1:2
    eps = epsv(j);
    yu0 = [eps 0];
    subplot(2,1,j)
    [t, yu] = ode45('pendulum_dydt',tspan,yu0);
    y1 = cos(t);
    y3 = (cos(t) - cos(3*t))/192 + t.*sin(t)/16;
    plot(t/(2*pi),yu(:,1),'b-',t/(2*pi),eps*y1,'r--', ...
         t/(2*pi),eps*y1 + eps^3*y3,'g--')
    xlabel('t/2\pi')
    ylabel('y(t) (radians)')
    title(['Nonlinear Pendulum, y(0) = \epsilon = ' num2str(eps)  ...
	   ', dy/dt(0) = 0'])
    legend('Exact','\epsilon y_1(t)', ...
          '\epsilon y_1(t) + \epsilon^3 y_3(t)')
  end
  hold off