%=======================================================================
% Chapter 10, code 1: a linear boundary value problem.
%
% Solve y'' + p(x)*y' + q(x)*y = r(x)
% with prescribed endpoint values y(a) = alpha and y(b) = beta.
%
% Method: centered finite differences on a uniform grid.
% Solve the resulting tridiagonal system with the Thomas algorithm.
%
% Example: y'' - y = 0, y(0) = y(1) = 1.
% Exact solution: cosh(x-0.5)/cosh(0.5).
%
% MATLAB alternative: assemble the interior system as a sparse matrix A,
% then solve with interior = A\rhs. The Thomas routine is retained here
% to show the elimination steps explicitly.
%
% To solve another problem:
%   1. Edit the interval, boundary values, and number of intervals.
%   2. Edit coefficients below the driver.
%   3. Disable check_exact or replace exact_solution.
%
% Run this file as a script. Local functions are placed at the end.
%=======================================================================
clear;

%-------------------- Problem and grid: edit here -----------------------
a = 0.0;
b = 1.0;
alpha = 1.0;                 % y(a)
beta = 1.0;                  % y(b)
n = 20;                     % n intervals, n+1 points, n-1 unknowns
check_exact = true;
make_plot = true;
%-----------------------------------------------------------------------
if n < 2 || n ~= floor(n) || b <= a
    error('Use an integer n >= 2 and b > a.');
end

h = (b-a)/n;
x = linspace(a,b,n+1)';
y = finite_difference(@coefficients,n,x,alpha,beta);

filename = 'ch10_finite_difference.dat';
fid = fopen(filename,'w');
if fid == -1
    error('Cannot open the data file.');
end

fprintf('Linear boundary value problem: finite differences\n');
fprintf('Grid intervals: %8d\n',n);
fprintf('Grid spacing:   %12.6f\n\n',h);

if check_exact
    exact = exact_solution(x,a,b,alpha,beta);
    abs_error = abs(y-exact);
    fprintf('%10s %16s %16s %12s\n', ...
            'x','numerical y','exact y','abs error');
    fprintf(fid,'# x numerical_y exact_y abs_error\n');
    for i = 1:n+1
        fprintf('%10.6f %16.8e %16.8e %12.4e\n', ...
                x(i),y(i),exact(i),abs_error(i));
        fprintf(fid,' %.16e %.16e %.16e %.16e\n', ...
                x(i),y(i),exact(i),abs_error(i));
    end
else
    fprintf('%10s %16s\n','x','numerical y');
    fprintf(fid,'# x numerical_y\n');
    for i = 1:n+1
        fprintf('%10.6f %16.8e\n',x(i),y(i));
        fprintf(fid,' %.16e %.16e\n',x(i),y(i));
    end
end
fclose(fid);

fprintf('\n');
if check_exact
    fprintf('Largest grid error: %12.4e\n',max(abs_error));
end
fprintf('Data file: %s\n',filename);

if make_plot
    figure;
    if check_exact
        subplot(2,1,1);
    end
    plot(x,y,'ko-','LineWidth',1.0,'MarkerSize',4);
    hold on;
    if check_exact
        x_plot = linspace(a,b,401)';
        plot(x_plot,exact_solution(x_plot,a,b,alpha,beta), ...
             'k--','LineWidth',1.2);
        legend('Finite differences','Exact solution','Location','best');
    end
    hold off;
    xlabel('x');
    ylabel('y(x)');
    title('Linear boundary value problem');
    grid on;

    if check_exact
        subplot(2,1,2);
        plot(x,abs_error,'ko-','LineWidth',1.0,'MarkerSize',4);
        xlabel('x');
        ylabel('Absolute error');
        grid on;
    end
end


function [p,q,r] = coefficients(x)
%=======================================================================
% Supply the coefficients of the original differential equation:
%   y'' + p(x)*y' + q(x)*y = r(x).
% Default example: y'' - y = 0.
%=======================================================================
p = 0.0;
q = -1.0;
r = 0.0;
end


function y = exact_solution(x,a,b,alpha,beta)
% Exact solution of y'' - y = 0 with the supplied endpoint values.
% Replace this function when changing the equation, or disable the check.
y = (alpha*sinh(b-x)+beta*sinh(x-a))/sinh(b-a);
end


function y = finite_difference(fcn,n,x,alpha,beta)
%=======================================================================
% Centered approximations at an interior point:
%   y'' = (y_left-2*y_center+y_right)/h^2
%   y'  = (y_right-y_left)/(2*h).
%
% Multiplying the differential equation by h^2 gives:
%   lower*y_left + diagonal*y_center + upper*y_right = h^2*r.
%
% MATLAB indices start at 1. Interior unknown i lies at x(i+1).
% Only the interior values are unknown. Move the known endpoint
% contributions to the right-hand side before solving.
%=======================================================================
h = (x(end)-x(1))/n;
m = n-1;
lower = zeros(m,1);
diagonal = zeros(m,1);
upper = zeros(m,1);
rhs = zeros(m,1);

for i = 1:m
    [p,q,r] = fcn(x(i+1));
    lower(i) = 1-h*p/2;
    diagonal(i) = -2+h*h*q;
    upper(i) = 1+h*p/2;
    rhs(i) = h*h*r;
end

% Both corrections apply when there is only one interior point.
rhs(1) = rhs(1)-lower(1)*alpha;
rhs(m) = rhs(m)-upper(m)*beta;
lower(1) = 0;
upper(m) = 0;

interior = thomas(lower,diagonal,upper,rhs);
y = [alpha; interior; beta];
end


function y = thomas(lower,diagonal,upper,rhs)
%=======================================================================
% Solve a tridiagonal system by forward elimination and backward
% substitution. The local copies of diagonal and rhs are overwritten.
%
% No pivoting is used. A near-zero pivot stops the calculation.
% Such a failure can require a solver with pivoting or can indicate
% a singular boundary value problem.
%=======================================================================
m = length(rhs);
y = zeros(m,1);

% Check each pivot before division.
for i = 1:m
    scale = max([abs(lower(i)),abs(diagonal(i)),abs(upper(i))]);
    if abs(diagonal(i)) <= eps*scale
        error('Thomas algorithm: zero or near-zero pivot.');
    end
    if i == m
        break;
    end
    multiplier = lower(i+1)/diagonal(i);
    diagonal(i+1) = diagonal(i+1)-multiplier*upper(i);
    rhs(i+1) = rhs(i+1)-multiplier*rhs(i);
end

y(m) = rhs(m)/diagonal(m);
for i = m-1:-1:1
    y(i) = (rhs(i)-upper(i)*y(i+1))/diagonal(i);
end
end