function m550bvp6

% The ODE:
%      y'' +  f(x) y = 0,  f = lambda^2 x^2 sin^2 x
% The BCs:
%      y(0) = y(pi/2) = 0  plus  y'(0) = 1 (normalization)
   

fun = []; lams=[]; lame=[];

% First we abandon the eigenvalue problem and shoot from x=0
% using y(0)=0 and y'(0)=1 as ICs, taking lambda as a varying
% parameter and monitoring y(pi/2)

for nlam=1:300;
    llam = nlam/200*15;                       % current lambda choice
    [t,y] = ode45(@odef,[0 pi/2],[0 1]);     % solve ODE by shooting
    fun = [fun y(end,1)]; lams=[lams llam];  % record results

% If the value of y(pi) switches sign, we're close to an eigenvalue;
% switch to boundary-value solver (BVP4C) and find it, using a quick
% finite difference to find a guess

    if nlam>2 & fun(nlam)*fun(nlam-1)<0
       dfun = (fun(nlam)-fun(nlam-1))/(lams(nlam)-lams(nlam-1));
       clear sol
       sol.parameters = lams(nlam)-fun(nlam)/dfun;
       sol.x=t; sol.y(1,:)=y(:,1)'; sol.y(2,:)=y(:,2)';
       sol = bvp4c(@ode,@bc,sol);
       subplot(211), plot(sol.x,sol.y(1,:)), hold on
       lame=[lame sol.parameters];
    end
end
xlabel('x')
ylabel('y(x)')
xlim([0 pi/2])
title('Eigenfunctions')
hold off

%  Number of eigenvalues detected
ns=1:length(lame);

% Naive WKB prediction
lamas=ns*pi;

% WKB prediction
lama=ns*pi-pi/6;

subplot(223)
plot(lams,fun,lame,lame*0,'*')
xlabel('\lambda')
ylabel('y(\pi/2)')
title('Shooting vs BVP4c')

subplot(224)
plot(ns,lamas,'s',ns,lama,'o',ns,lame,'*')
lame,lamas,lama
ylabel('\lambda')
xlabel('n')
title('WKB vs BVP4c')

function f = funk(x,lam)
  f = lam^2*x^2*sin(x)^2;
end

function dydx = odef(x,y)         % ODE for shooting
  dydx = [              y(2)
  	       -funk(x,llam)*y(1) ];
end 

function dydx = ode(x,y,lambda)   % ODE for eigenvalue problem
  dydx = [              y(2)
	       -funk(x,lambda)*y(1) ];
end
   
function res = bc(ya,yb,lambda)   % BCs for eigenvalue problem
  res = [  ya(1)
           yb(1)
          ya(2)-1 ];
end 

end
