%
%  code16.m -- duffing oscillator (22 oct 2001 - djm)
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global A B period
A = 1.0;  B = (0.1);  tol = 1e-8;  period = 2*pi;

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

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

Xode = 'code16a';  IVs = [U0 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

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

[t y] = feval(method,Xode,[0 T],IVs,options);
amp = 2*y(end,3)/period;

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%  solution plot

clf;subplot(2,1,1)

plot(t,y(:,1),'kx')
hold on
plot(t,y(:,1),'b')
plot(t,y(:,2),'kx')
plot(t,y(:,2),'r')

title('\bf numerical duffing oscillation (\epsilon = 0.1)')
xlabel('\bf t-axis');  ylabel('\bf y(t) in blue, y\prime(t) in red');

amax = max(max(abs(y)));  Tend = max(t);  
axis([0 Tend -1.1*amax 1.1*amax])

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%  asymptotic solution plot

y_ex = y(:,1);
y_as = amp*cos(2*pi*t/period) - (B/32)*(amp^3)*cos(3*2*pi*t/period);
y_as = y_as + (B^2/32)*(amp^5)*cos(3*2*pi*t/period)/16;
y_as = y_as + (B^2/32)*(amp^5)*cos(5*2*pi*t/period)/96;
y_as = y_as - (B^2*27/32/32)*(amp^5)*cos(3*2*pi*t/period);

subplot(2,1,2)
plot(t,y_ex-y_as,'b')
hold on
%plot(t,y(:,1),'k')
hold off

title('\bf asymptotic error (\epsilon = 0.1)')
xlabel('\bf t-axis');  ylabel('\bf y_{num}(t)-y_{asy}(t)');

%  first integral

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