%=======================================================================
% Chapter 10, code 4: eigenvalues by Numerov shooting and bisection.
% Solve y'' = g(x,lambda)*y, with y(a) = y(b) = 0.
% Scan lambda for endpoint sign changes, then refine each bracket.
%
% Example: a uniform vibrating string, y'' = -lambda*y on [0,1].
% Exact eigenvalues: lambda_m = (m*pi/(b-a))^2, m = 1,2,... .
% Normalize each mode so that the integral of y^2 over [a,b] is 1.
%
% MATLAB alternatives: fzero can refine a sign-changing shooting bracket.
% A finite-difference matrix formulation can use eig or sparse eigs;
% this changes the discretization from Numerov shooting.
%
% To solve another problem:
%   1. Edit the interval, grid, number of modes, and search controls.
%   2. Edit coefficient below the driver.
%   3. Disable check_exact or replace the analytic comparisons below.
% This version assumes zero-value boundaries at BOTH ends and no y' term.
%=======================================================================
clear;

%-------------------- Problem and search: edit here ----------------------
a = 0.0;
b = 1.0;
n = 200;                     % Number of grid intervals
modes = 3;
lambda_min = 0.0;
lambda_max = 100.0;
scan_step = 1.0;
tolerance = 1.0e-10;          % Maximum width of the final lambda bracket
check_exact = true;
make_plot = true;
% Increase lambda_max if the range contains too few modes.
% Make scan_step small enough to separate consecutive roots.
%------------------------------------------------------------------------
if n < 2 || n ~= floor(n) || modes < 1 || modes ~= floor(modes) || b <= a
    error('Use integer n >= 2, integer modes >= 1, and b > a.');
end
if lambda_max <= lambda_min || scan_step <= 0 || tolerance <= 0
    error('Use an increasing search interval, positive scan step and tolerance.');
end
h = (b-a)/n;
x = linspace(a,b,n+1)';

fprintf('Eigenvalues: Numerov shooting and bisection\n');
fprintf('Grid intervals:    %8d\n',n);
fprintf('Grid spacing:      %12.6f\n',h);
fprintf('Bracket tolerance: %12.4e\n',tolerance);
[values,shapes,residuals,iterations] = find_modes(@coefficient,n,x,modes, ...
    lambda_min,lambda_max,scan_step,tolerance);

fprintf('\n');
if check_exact
    fprintf('%6s %18s %18s %12s %10s\n', ...
        'mode','numerical lambda','exact lambda','abs error','bisections');
else
    fprintf('%6s %18s %10s\n','mode','numerical lambda','bisections');
end
for m = 1:modes
    if check_exact
        exact = (m*pi/(b-a))^2;
        fprintf('%6d %18.10f %18.10f %12.4e %10d\n', ...
            m,values(m),exact,abs(values(m)-exact),iterations(m));
    else
        fprintf('%6d %18.10f %10d\n',m,values(m),iterations(m));
    end
end

% One file for all modes. Each column of shapes contains one eigenfunction.
filename = 'ch10_numerov_eigenvalues.dat';
fid = fopen(filename,'w');
if fid == -1
    error('Cannot open the data file.');
end
fprintf(fid,'# Numerov eigenfunctions: x y_1 y_2 ...\n');
fprintf(fid,'# Each mode has integral(y^2 dx) = 1 (trapezoidal rule).\n');
for m = 1:modes
    fprintf(fid,'# mode %6d lambda = %.16e\n',m,values(m));
end
fprintf('\nEigenfunctions: columns are x, then the normalized modes.\n');
shape_errors = zeros(modes,1);
for i = 1:n+1
    fprintf('%10.6f',x(i));
    for m = 1:modes
        fprintf(' %16.8e',shapes(i,m));
        if check_exact
            exact_shape = sqrt(2/(b-a))*sin(m*pi*(x(i)-a)/(b-a));
            shape_errors(m) = max(shape_errors(m),abs(shapes(i,m)-exact_shape));
        end
    end
    fprintf('\n');
    fprintf(fid,' %.16e',[x(i),shapes(i,:)]);
    fprintf(fid,'\n');
end
fclose(fid);
fprintf('\n%6s %18s\n','mode','endpoint residual');
for m = 1:modes
    fprintf('%6d %18.8e\n',m,residuals(m));
end
if check_exact
    fprintf('%6s %20s\n','mode','largest shape error');
    for m = 1:modes
        fprintf('%6d %20.8e\n',m,shape_errors(m));
    end
end
fprintf('Data file: %s\n',filename);

if make_plot
    figure;
    x_plot = linspace(a,b,801)';
    for m = 1:modes
        subplot(modes,1,m);
        plot(x,shapes(:,m),'k-','LineWidth',1.2);
        hold on;
        if check_exact
            exact_shape = sqrt(2/(b-a))*sin(m*pi*(x_plot-a)/(b-a));
            plot(x_plot,exact_shape,'k--','LineWidth',1.0);
            legend('Numerov','Exact solution','Location','best');
        end
        hold off;
        xlabel('x');
        ylabel('y_m(x)');
        title(sprintf('Mode %d: lambda = %.8f',m,values(m)));
        grid on;
    end
end


function g = coefficient(x,lambda)
% Supply g(x,lambda) in y'' = g(x,lambda)*y.
% Uniform string: g = -lambda.
% For variable dimensionless density rho(x), use g = -lambda*rho(x).
% Numerov requires a linear homogeneous equation with no y' term.
g = -lambda;
end


function [values,shapes,residuals,iterations] = find_modes(fcn,n,x,modes, ...
    lambda_min,lambda_max,scan_step,tolerance)
% The shooting mismatch is y(b;lambda). Eigenvalues make it zero.
% A sign scan can miss two or more roots inside one interval. Reduce
% scan_step and repeat when adapting this code to another equation.
% Analytic mode numbers assume the scan finds the lowest modes in order.
values = zeros(modes,1);
shapes = zeros(n+1,modes);
residuals = zeros(modes,1);
iterations = zeros(modes,1);
found = 0;
left = lambda_min;
[trial,f_left] = numerov(fcn,n,x,left);
% Handle an eigenvalue exactly at the initial scan endpoint once.
if f_left == 0
    found = 1;
    values(found) = left;
    trial = normalize_mode(n,x,trial);
    shapes(:,found) = trial;
    residuals(found) = trial(end);
end

while found < modes && left < lambda_max
    right = min(left+scan_step,lambda_max);
    if right <= left
        error('Scan step cannot advance lambda.');
    end
    [~,f_right] = numerov(fcn,n,x,right);
    bracketed = (f_left < 0 && f_right > 0) || (f_left > 0 && f_right < 0);
    if f_right == 0 || bracketed
        if f_right == 0
            root = right;
            count = 0;
        else
            [root,count] = bisect(fcn,n,x,left,right,f_left,tolerance);
        end
        found = found+1;
        values(found) = root;
        iterations(found) = count;
        [trial,~] = numerov(fcn,n,x,root);
        trial = normalize_mode(n,x,trial);
        shapes(:,found) = trial;
        % Keep the normalized endpoint residual; do not force it to zero.
        residuals(found) = trial(end);
    end
    left = right;
    f_left = f_right;
end
if found < modes
    error('Too few modes found: enlarge the range or reduce scan_step.');
end
end


function [root,iterations] = bisect(fcn,n,x,left,right,f_left,tolerance)
% Stop on lambda-bracket width, not on mismatch magnitude, which depends
% on the arbitrary shooting amplitude. Preserve the scan endpoints.
lo = left;
hi = right;
flo = f_left;
iterations = 0;
max_iterations = 100;
while hi-lo > tolerance
    if iterations >= max_iterations
        error('Bisection iteration limit reached.');
    end
    mid = lo+(hi-lo)/2;
    if mid <= lo || mid >= hi
        error('Lambda bracket cannot shrink; check tolerance.');
    end
    [~,fmid] = numerov(fcn,n,x,mid);
    iterations = iterations+1;
    if fmid == 0
        root = mid;
        return;
    end
    if (flo < 0 && fmid > 0) || (flo > 0 && fmid < 0)
        hi = mid;
    else
        lo = mid;
        flo = fmid;
    end
end
root = lo+(hi-lo)/2;
end


function [y,mismatch] = numerov(fcn,n,x,lambda)
%=======================================================================
% Numerov for y'' = g*y, with c = h^2/12:
%   (1-c*g_next)*y_next =
%       2*(1+5*c*g_center)*y_center - (1-c*g_previous)*y_previous.
%
% MATLAB indices start at 1. y(1) = 0 and y(2) = h fix an arbitrary
% nonzero amplitude. No Taylor approximation of a derivative is needed:
% the homogeneous recurrence makes eigenvalues independent of amplitude.
% Normalize the completed mode after finding the root.
% This starting choice applies to a zero-value left boundary only.
%=======================================================================
h = (x(end)-x(1))/n;
c = h*h/12;
g = zeros(n+1,1);
y = zeros(n+1,1);
for i = 1:n+1
    g(i) = fcn(x(i),lambda);
    if ~isfinite(g(i))
        error('Nonfinite Numerov coefficient.');
    end
end
y(1) = 0;
y(2) = h;
for i = 2:n
    denominator = 1-c*g(i+1);
    if abs(denominator) <= eps*max(1,abs(c*g(i+1)))
        error('Near-zero Numerov denominator: refine the grid.');
    end
    y(i+1) = (2*(1+5*c*g(i))*y(i)-(1-c*g(i-1))*y(i-1))/denominator;
    if ~isfinite(y(i+1))
        error('Nonfinite Numerov solution: check equation, interval, and grid.');
    end
end
mismatch = y(end);
end


function y = normalize_mode(n,x,y)
% Trapezoidal approximation to integral y^2 dx; keep the initial sign.
h = (x(end)-x(1))/n;
norm_value = sqrt(h*(0.5*y(1)^2+sum(y(2:n).^2)+0.5*y(end)^2));
if norm_value <= 0 || ~isfinite(norm_value)
    error('Cannot normalize the eigenfunction.');
end
y = y/norm_value;
end