% author: Akbar Alvi Rahman
% Lecture 1 Worked Example 1 (Section 10.1)

clear all; close all; clc

% define physical constants
L = 1;
S = 0;
conductivity = 400;
Ta = 300;
Tb = 320;

Sp = -100;
Sc = 5000;

% define numerical parameters
n = 21;

% grid generation
x0 = linspace(0, L, 21);
dx = L/(n-1);
Dx = dx;
DxB = Dx/2;

% create matrix
A = zeros(n, n);
B = zeros(n, 1);

A(1, 1) = 1; % Dirichlet boundary
B(1) = Ta;

A(n, n) = 1; % Dirichlet boundary
B(n) = Tb;

for i=2:n-1
    A(i, i-1) = -conductivity/dx;
    A(i, i) = 2*conductivity/dx-Sp*Dx;
    A(i, i+1) = -conductivity/dx;

    B(i) = Sc*Dx;
end

T = A\B; % matrix backslash

figure('color', 'w'); 
plot(x0, T, 'rs');
grid on;
xlabel('x [m]'); ylabel('T [K]');
hold on;

% comparison with theoretical solution
Tteo = Ta + (Tb-Ta)*x0/L;
plot(x0, Tteo, 'k-');
legend('backslash', 'Exact');
error = mean(abs(T-Tteo')) % Tteo transposed by ' operator

% check of residuals
residual = sum(abs(B-A*T))/sum(abs(diag(A).*T)) % defined as |B-A*T|/|diag(A).T|

% check of energy balance
num=0; den = 0;

for i=2:n-1
    num = num+abs(-conductivity*(T(i)-T(i-1))/dx+conductivity*(T(i+1)-T(i))/dx); % sum(|q_w-q_e|)
    den = den+abs(conductivity*(T(i)-T(i-1))/dx); % sum(|q_w|)
end

unbalance = num/den*100