%
%  generalized WKB -- code20.m (23 oct 2001 - djm)
%
%- - - - - - - - - - - - - - - - - - - - - - - - 

%  coordinate matrices

du = 0.02;  
u  =    du/2:du:2+du/2;
up = -1.5-du/2:du:1.5+du/2;  [u1 u2] = meshgrid(u,up);

%  first integral w/o damping

global kk;  kk = 1;

integ = (u2.^2)/2 + (kk^2)*(u1.^2)/2 + (1/2)./(u1.^2);
c0   = 0.5:0.1:0.9;
cont = (kk^2)*(c0.^2)/2 + (1/2)./(c0.^2);

%  contourplot

figure(2);  clf;  
figure(1);  clf;  
contour(u,up,integ,cont);  hold on;
title('\bf amplitude phase plane');
xlabel('\bf A(t)');  ylabel('\bf A\prime(t)')
axis([0 2 -2 2]); axis square

while(1)

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

  Xode = 'code20a';  T = 100;  method = 'ode45';  tol = 1e-6;

%  input IVs

  figure(1);  clear IVs;  

  IVs = input('initial A (0<A<1):  ');

  if(isempty(IVs))
    break
  else
    IVs = [IVs 0 0];
    plot(IVs(1),IVs(2),'kx')
  end

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

%  ODE solve & plot

  [t y] = feval(method,Xode,[0 T],IVs,options);
  plot(y(:,1),y(:,2),'k.','markersize',10)
  disp(['period:  ' num2str(y(end,3))])

  figure(2);  
  subplot(2,1,1);  hold on
  plot(t,cos(y(:,3)).*y(:,1))
  title('\bf original sinusoids')
  xlabel('\bf t-axis');  ylabel('\bf y(t)')
  
  subplot(2,1,2);  
  plot(t,y(:,1),'b');  hold on
  plot(t,y(:,2),'r');  hold off
  title('\bf amplitude solutions')
  xlabel('\bf t-axis');  ylabel('\bf A(t)=blue, A\prime(t)=red')
end