%
%  code15.m -- duffing oscillator
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global A B period
A = 1.0;  B = 0.1;  tol = 1e-12;  period = 2*pi;

%  ODE file, max time, method & options
Xode = 'code15a';  T = 10;  method = 'ode45';

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

%  parametric ODE solve
U0 = 2.^(-5:0.5:1);  per = [];  amp = [];

for j=1:size(U0,2)
  IVs = [U0(j) 0 0];

%  period solve
  [t y] = feval(method,Xode,[0 T],IVs,options);

%  amplitude integral
  period = 4*t(end);
  [t y] = feval(method,Xode,[0 T],IVs,options);

%  first integral check
  check = (y(:,2).^2)/2 + A*(y(:,1).^2)/2 - B*(y(:,1).^4)/4 ;
%  disp(['1st-integral variation = ' num2str(max(check)-min(check))])

  amp = [amp ; 8*y(end,3)/period];
  per = [per ; period-2*pi];
  
  data = [U0(j) amp(j) per(j)]
end

per2 = 2*pi./sqrt(1 - (3/4)*B*amp.^2) - (2*pi);
per4 = 2*pi./sqrt(1 - (3/4)*B*amp.^2 + (3/128)*(B^2)*amp.^4) - (2*pi);

max0 = amp;
max2 = max0 - (1/32)*B*(amp.^3);
max4 = max2 + (3/32)*(B^2).*(amp.^5)*((1/16)+(1/96)-(9/32)) ;

clf
subplot(2,1,1)
loglog(amp,abs(per),'b+')
hold on
loglog(amp,abs(per2-per),'kx')
loglog(amp,abs(per4-per),'ro')
plot([0.1 1],[1e-2 1e0],'b')
plot([0.1 1],[1e-6 1e-2],'k')
plot([0.1 1],[1e-9 1e-3],'r')
hold off
axis([2^(-5) 2.5 1e-14 1e1])
title('\bf period error (slopes = 2,4,6)')
xlabel('\bf amplitude of cos(\tau)-mode')
ylabel('\bf |P_{num}-P_{asy}|')
text(0.3,1e-12,'\bf \epsilon = 0.1')

subplot(2,1,2)
loglog(amp,abs(max0-U0'),'b+')
hold on
loglog(amp,abs(max2-U0'),'kx')
loglog(amp,abs(max4-U0'),'ro')
plot([0.1 1],[1e-5 1e-2],'b')
plot([0.1 1],[1e-8 1e-3],'k')
plot([0.1 1],[1e-11 1e-4],'r')
hold off
axis([2^(-5) 2.5 1e-14 1e-0])
title('\bf maximum error (slopes = 3,5,7)')
xlabel('\bf amplitude of cos(\tau)-mode')
ylabel('\bf |y_{num}(0,\epsilon)-y_{asy}(0,\epsilon)|')
text(0.3,1e-12,'\bf \epsilon = 0.1')