%
%  generalized WKB -- code18.m (22 oct 2001 - djm)
%
%- - - - - - - - - - - - - - - - - - - - - - - - 

%  coordinate matrices

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

%  first integral w/o damping

global kk;  kk = 1;

integ = (u2.^2)./(u1.^3) + 4*u1 + 4*(kk^2)./u1;

%  contourplot

figure(2);  clf;  
figure(1);  clf;  
contour(u,up,integ,[8.2 9.2 10.2 12.2]);  hold on;
contour(u,up,integ,-[8.2 9.2 10.2 12.2]);  cbar
title('\bf frequency phase plane');
xlabel('\omega(t)');  ylabel('\omega\prime(t)')

while(1)

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

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

%  input IVs

  figure(1);  clear IVs;  

  IVs = input('initial value (0<omega<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','on','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))./sqrt(abs(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 frequency solutions')
  xlabel('\bf t-axis');  ylabel('\omega(t)=blue, \omega\prime(t)=red')
end