%
%  code17.m -- slowly-varying oscillator
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global ee
ee = 0.25;  tol = 1e-8;  T0 = 0;

for j=1:8
  j
  
%  ODE file, IVs, final time, method & options
  Xode = 'code17a';  IVs = [1 0];  T = (2*pi)/(ee^2);  method = 'ode45';

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

%  ODE solve 
  [t y] = feval(method,Xode,T0+[0 T],IVs,options);
  IVs = [y(end,:)];

%  phase integration
  Xode = 'code17b';  IVph = [0 0];  T = (2*pi)/(ee^2);  method = 'ode45';

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

  [tph ph] = feval(method,'code17b',T0+[0 T],[0 0],options);
  IVph = [ph(end,:)];  T0 = t(end);

%  amplitude

  k0 = 2 + cos(ee*tph);
  k1 = -ee*sin(ee*tph);
  k2 = -(ee^2)*cos(ee*tph);

  s = k0 + (3/8)*(k1.^2)./(k0.^3) - (1/4)*(k2)./(k0.^2);
end

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

clf;subplot(2,1,1)
plot(t,y(:,1),'b')
hold on
plot(tph,sqrt(s(1)./s),'r')
plot(tph,-sqrt(s(1)./s),'r')
plot(tph,sqrt(k0(1)./k0).*cos(ph(:,1)),'g')
plot(tph,sqrt(s(1)./s).*cos(ph(:,2)),'k')

title('\bf long-time solution (y_{num}=blue, y_{leading}=green, y_{1st corr}=black/red)')
xlabel('\bf t-axis');  ylabel('\bf y-axis');

amax = max(abs(y(:,1)));
axis([T0-T T0 -1.1*amax 1.1*amax])

subplot(2,1,2)
plot(t,y(:,1),'b')
hold on
plot(tph,sqrt(s(1)./s),'r')
plot(tph,-sqrt(s(1)./s),'r')
plot(tph,sqrt(k0(1)./k0).*cos(ph(:,1)),'gx')
plot(tph,sqrt(s(1)./s).*cos(ph(:,2)),'kx')

title('\bf long-time solution (phase slip)')
xlabel('\bf t-axis');  ylabel('\bf y-axis');

amax = max(abs(y(:,1)));
axis([T0-ee*T T0 -1.1*amax 1.1*amax])