%
%  finder.m -- (y lee:  19 mar 2004)
%    - sturm-liouville eigenvalue code
%

clear

N = 50;  fig_N = 10;

AA = 1; BB = 500;
lambda = linspace(sqrt(AA), sqrt(BB), N).^2;

a = 1; b = 2;
ua = 0;
ode_rel_tol = 1e-7; ode_abs_tol = 1e-10;
fzero_opt = optimset('TolFun',1e-10);

for j = 1:N
	ub(j) = u_ivp(lambda(j), a, b, ua, ode_rel_tol, ode_abs_tol,2);
end

%  plot u(b,lambda)
figure(fig_N);
hold off
plot(lambda, ub,'k.')
axis( [AA BB min(ub) max(ub) ] )
grid on;  hold on

title(['\bf u(b;\lambda)'])
xlabel(['\bf \lambda'])
ylabel(['\bf zeros are eigenvalues'])

% Let's search!  look for sign changes
pr = ub(1:N-1).*ub(2:N);
lambda_zero = lambda( pr < 0);
plot( lambda_zero, 0, 'o');

num_lambda_to_search = length(lambda_zero);

for k = 1:num_lambda_to_search
	my_lambda(k) = fzero(@u_ivp, lambda_zero(k), fzero_opt, a, b, ua, ...
			ode_rel_tol, ode_abs_tol,2);
end
plot( my_lambda, 0, '*');

my_lambda

%  eigenfunction plot
ode_params = [a, b, ua]
ode_tols = [ode_rel_tol, ode_abs_tol]

for k = 1: num_lambda_to_search
	[ub, u_k, x_k] = u_ivp(my_lambda(k), a, b, ua, ...
                       ode_rel_tol, ode_abs_tol,200);
	u_all(:,k) = u_k;
end

figure(fig_N+1);
plot(x_k, u_all)

title(['\bf first ' num2str(num_lambda_to_search) ' eigenfunctions'])
xlabel(['\bf x-axis'])
ylabel(['\bf u(x)'])