%
%  code22.m -- limit cycles (31 oct 2001 - djm)
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global epsil
epsil = 0.2;  tol = 1e-6;

disp(' ')
U0 = input('initial derivative:  ');

%  ODE file, IVs, final time, method & options

Xode = 'code22a';  IVs = [0 U0];  T0 = 0.0;  T = 10;  method = 'ode45';

options = odeset('reltol',tol,'abstol',tol,'refine',1, ...
		'maxstep',T/20,'stats','off','events','on');

%  ODE solve with period & amplitude determination

figure(1);  axis equal;  hold on
figure(2);  hold on

for N = 1:8
  [t y] = feval(method,Xode,[T0 T0+T],IVs,options);
  disp(['period = ' num2str(t(end)-T0)]);

%  phase plane

  figure(1)
  plot(y(:,1),y(:,2),'b')
  plot(y(:,1),y(:,2),'k.')

  title(['\bf phase plane (\epsilon = ' num2str(epsil) ')'])
  xlabel('\bf y(t)');  ylabel('\bf y\prime(t)');  pause(0)

  IVs = [y(end,1),y(end,2)];  T0 = t(end);

%  plot
  figure(2)
  plot(t,sqrt(y(:,1).^2 + y(:,2).^2),'b')

  title(['\bf oscillations (\epsilon = ' num2str(epsil) ')'])
  xlabel('\bf t-axis');  ylabel('\bf sqrt(y(t)^2 + y\prime(t)^2)');
end