%{
  Solving ODE Boundary value problem 
  y'' + P(x)y' +Q(x)y = F(x)
  Boundary conditions: chose by key value
    1. Dirichlet a,  y(a), b,  y(b)
    2. mixed as  a,  y(a), b, y'(b)
    3. mixed as  a, y'(a), b,  y(b)
  Method: Finite Dimmerece method (equilibrium)
  AG: 2 April 2022
  ToDo: add 4th boundary condition a, y'(a), b, y'(b)
%}

% 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

% n grid size
N = 41;

% boundary conditions (Dirichlet): a, y(a), b, y(b)
% key - type of boundary conditions
%   1 - Dirichlet
%   2 - mixed as y(a), y'(b)
%   3 - mixed as y'(a), y(b)
key = 3;
a  =  0.0;
ya =  0.0;
b  =  2.0;
yb =  0.0;
% functions in y'' + P(x)y' + Q(x)y = F(x)
P = @Px;
Q = @Qx;
F = @Fx;

% ==== solver

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

% ====

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

%=== 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;
Px = 2.0;
end

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

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

%{
Solving linear boundary value problem with various 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
key - type of boundary conditions
  1 - Dirichlet
  2 - mixed as y(a), y'(b)
  3 - mixed as y'(a), y(b)
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)

AG
%}

function [X,Y] = ODEbvE2(P, Q, F, key, 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;           % ai in my notes
    c(k,2) = (1.0)*h*h*Q(x(k)) - 2.0;       % bi in my notes
    c(k,3) = 1.0 + h*P(x(k))/2.0;           % ci in my notes
    d(k)   = h*h*F(x(k));                   % 
end

switch key
    case 1
        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;
    case 2
        d(1)   = d(1) - c(1,1)*ya;
        c(1,1) = 0.0;
        d(n)   = d(n) - 2.0*h*c(n,3)*yb;
        c(n,1) = 2.0;
        c(n,3) = 0.0;
    case 3
        d(1)   = d(1) - 2.0*h*c(1,1)*ya;
        c(1,1) = 0.0;
        c(1,3) = 2.0;
        d(n)   = d(n) - c(n,3)*yb;
        c(n,3) = 0.0;
end

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

switch key
    case 1
        X(1)  = a;
        X(N)  = b;
        Y(1)  = ya;
        Y(N)  = yb; 
    case 2
        X(1)  = a;
        X(N)  = b;
        Y(1)  = ya;
        Y(N)  = y(n) + h*yb; 
    case 3
        X(1)  = a;
        X(N)  = b;
        Y(1)  = y(1) - 1.0*h*ya;
        Y(N)  = yb;         
end

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