function [ lambda_num, lambda_WKB, sol ] = WKBex( n, solver )
%WKBex  Use bvp4/5c() to compute the nth eigval and assoc'd eigfun of the EVP
% 
%         phi'' + lambda*(1+x)*phi = 0,  0 < x < 1,  phi(0) = phi(1) = 0
% 
%       Example pp. 207-8 [Haberman, "Elem Appl PDEs w/ Fourier Series and
%       BVPs", 3rd ed ] and compare against WKB approx.
% 
% >> [ lambda_num, lambda_WKB, sol ] = WKBex( n, solver )
% 
%   in:       n = index of eigen mode to approx (default = 1)
%        solver = which BVP solver to use: @bvp4c (old) or @bvp5c (new, default)
% 
%  out:  lambda_num = numerical approx eigval
%        lambda_WKB = WKB approx eigval
%               sol = struct containing numerical approx sol'n from solver:
%                  .x  = x(1:npts) grid (generated/adapted in solver)
%                  .y  = y(1:npts) approx eigfun ( phi, phi' ) on grid
%                  .yp = yp(1:npts) approx deriv's of y vect on grid
%                  .parameters = lambda = numer approx eigval
% 
% also output info to terminal window + plot

% ~/matlab/courses/AppliedMath/WKBex.m
% 
% C.G. (10-03)
% C.G. (11-09) changed to use either bvp4c() to bvp5c(), plus minor edits

% set params/defaults

if ( nargin < 1 || isempty(n) )
    n = 1 ;                            % default eigenmode = 1
end

if ( nargin < 2 || isempty(solver) )
    solver = @bvp5c ;                  % default solver = @bvp5c
end

% build initial guess

x = linspace( 0, 1, 4*n+1 ) ;          % 8 grid cells per period init grid

y = zeros( 2, length(x) ) ;            % allocate space for init sol'n guess

% common argument that goes into the trig functions

arg = n * pi * ( (1+x).^(1.5)-1 ) / ( sqrt(8) - 1 ) ;

y(1,:) = (1+x).^(-.25) .* sin(arg) ;   % init guess for phi (from WKB approx)

% init guess for phi' (from d/dx of WKB approx)

y(2,:) = 1.50 * n * pi * (1+x).^(.25) / ( sqrt(8)-1 ) .* cos(arg) ...
       -  .25 * (1+x).^(-1.25) .* sin(arg) ;

% init guess for eigval (from WKB approx)

lambda_WKB = 2.25 * n^2 * pi^2 / ( sqrt(8) - 1 )^2 ;

% put into init-guess struct

solinit = struct( 'x', x, 'y', y, 'parameters', lambda_WKB ) ;

% set options
% 
% note--with 'bvp4c', RelTol = 1e-11 and AbsTol = 1e-9 was about as far as
% I could push it with my examples (eigen modes 1-7) ...
% can go farther with 'bvp5c'

options = bvpset( ...
    'RelTol', 1e-11, ...               % {1e-3}
    'AbsTol', 1e-09, ...               % {1e-6}
    'Nmax',   10000, ...               % {floor(1000/n)}, wh n=ode syst size
    'Stats',  'on'  ) ;                % on | {off}

% wrapper function to pass additional argument to @bcfun ...
% by default, @odefun and @bcfun are called via odefun( x, y, lambda ) and
% bcfun( y0, y1, lambda ) --- bcfun doesn't need lambda but does need 'n'

my_bcfun = @( y0, y1, lambda ) bcfun( y0, y1, n ) ;

% call the solver

sol = solver( @odefun, my_bcfun, solinit, options ) ;

lambda_num = sol.parameters ;

% plot numerical vs WKB approx eigfuns

xi = linspace( 0, 1, 32*n+1 ) ;        % finer grid for plotting

yi = deval( sol, xi ) ;                % interpolate numerical sol'n onto xi

% WKB approx evaluated on fine grid

yWKB = (1+xi).^(-.25) .* sin( n*pi*( (1+xi).^(1.5)-1 ) / ( sqrt(8)-1 ) ) ;

clf                                    % clear current figure window

% plot numerical approx as blue solid, WKB approx as red dashed

plot( xi, yi(1,:), 'b-', xi, yWKB, 'r--' )

set( gca, 'FontSize', 12 )

legend( 'numerical', 'WKB' )           % add a legend

hold on

plot( [ 0 1 ], [ 0 0 ], 'k-' )         % add a horizontal axis (black solid)

axis( [ 0 1 -1 1 ] )

box on

xlabel( 'x' )
ylabel( '\phi' )

title_str = [ 'n = ', num2str(n), ...
        ',   \lambda_{num} = ', num2str(lambda_num,10), ...
        ',   \lambda_{WKB} = ', num2str(lambda_WKB,10) ] ;

title( title_str, 'FontSize', 14 ) ;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dydx = odefun( x, y, lambda )

dydx = [               y(2)
         -lambda*(1+x)*y(1) ] ;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function res = bcfun( y0, y1, n )

res = [ y0(1)
        y1(1)
        y0(2) - 1.5*n*pi/(sqrt(8)-1) ] ;

end