%{
  Solving ODE Boundary value problem for
  y'' + P(x)y' +Q(x)y = F(x)
  Boundary conditions - Dirichlet: a, y(a), b, y(b)
  Method: Finite Difference method (equilibrium)
  AG: April 2022
%}

% Matlab specific 
clear  % Clears command window
clc    % Clears command history
clf    % Removes anything in the figure window before simulation. 
tic    % set timer for the program

% INPUT ==========
N = 52;  % grid size 

% boundary conditions (Dirichlet): a, y(a), b, y(b)
a  =  0.0;
ya =  1.0;
b  =  1.0;
yb =  0.0;

% functions in y'' + P(x)y' + Q(x)y = F(x)
P = @Px;
Q = @Qx;
F = @Fx;
% end INPUT ==========

% ==== ODE solver

[X,Y] = ODEbvE(P, Q, F, a, ya, b, yb, N);

% ====

% plot solutions
figure (1)
plot(X,Y,'r-o')
xlabel('x')
ylabel('y(x)')
title('Boundary Value Problem');

%=== print solutions (if needed)
% fprintf( 'Solutions: \n       x            y \n')
% for k = 1:N
%     fprintf('%12.6f %12.6f \n', X(k),Y(k))
% end

%print time of calculations
timeE = ceil(toc);
fprintf('\n Elapsed time = %6i sec \n',timeE)


function Px = Px(x)
Px = sqrt(x)+5.0;
end

function Qx = Qx(x)
Qx = x-1.0;
end

function Fx = Fx(x)
Fx = cos(x) * exp(-x/4.0)  + 1.0;
end

%{
Solving linear boundary value problem with Dirichlet boundary conditions
y'' + P(x)y' +Q(x)y = F(x)

METHOD: Equilibrium + Thomas algorithm for solving 3-diagonal system

CALLS: external functions: Thomas.m

INPUT:
P(x), Q(x), F(x) - external functions
a and y(a), b and y(b) - boundary conditions
N - total number of points in the grid, including two boundary points

OUTPUT: 
Y(X) - solution as arrays Y and X (size N)

Last revision: AG April 2022
%}
function [X,Y] = ODEbvE(P, Q, F, a, ya, b, yb, N)
n = N-2;
h = (b-a)/(n+1);

c = zeros(n,3);
d = zeros(n,1);
x = zeros(n,1);

for k = 1:n
    x(k) = a + k*h;
    c(k,1) = 1.0 - h*P(x(k))/2.0;
    c(k,2) = (1.0)*h*h*Q(x(k)) - 2.0;
    c(k,3) = 1.0 + h*P(x(k))/2.0;
    d(k)   = h*h*F(x(k));
end

d(1)   = d(1) - c(1,1)*ya;
c(1,1) = 0.0;
d(n)   = d(n) - c(n,3)*yb;
c(n,3) = 0.0;

[y] = Thomas(c,d,n);

X(1)  = a;
X(N)  = b;
Y(1)  = ya;
Y(N)  = yb;
for k = 1:n
    X(k+1) = x(k);
    Y(k+1) = y(k);
end

end

function [x] = Thomas(c,b,n)
% !============================================================
% ! Solutions to a system of tridiagonal linear equations c*x=b
% ! Method: the Thomas method
% ! Alex G. March 2021
% !-----------------------------------------------------------
% ! input ...
% ! c(n,3) - array of coefficients for matrix where c(n,2) - diagonal
% ! b(n)   - vector of the right hand coefficients b
% ! n      - number of equations
% ! output ...
% ! x(n)   - solutions
% ! comments ...
% ! the original arrays c(n,3) and b(n) will be destroyed 
% ! during the calculation
% !===========================================================

%step 1: forward elimination
for k = 2:n
  coeff=c(k,1)/c(k-1,2);
  c(k,2)=c(k,2)-coeff*c(k-1,3);
  b(k)=b(k)-coeff*b(k-1);
end

%step 2: back substitution
x(n) = b(n)/c(n,2);
for k = n-1:-1:1
    x(k) = (b(k)- c(k,3)*x(k+1))/c(k,2);
end

end % Thomas