%
% Finite Difference method example of solving boundary value ODE
%
% d^2 T
% ----- + c (Tinf - T) = 0
% dx^2
%
% which models the temperature distribution in a rod
% define constants
Ta = 40; % left side temperature
Tb = 200; % right side temperature
Tinf = 20; % ambient temperature
c = 0.01; % constant
L = 10; % length of bar
ndivs = 6; % 'Input number of divisions;
nunknowns = ndivs - 1;
deltax = L/ndivs;
A = -(2 + deltax^2*c);
B = -deltax^2*c*Tinf;
matrix = zeros(nunknowns); % initialize
matrix(1,1) = A; % first row
matrix(1,2) = 1;
rhs(1) = B - Ta;
for i = 2:nunknowns - 1 % intermediate rows
matrix(i,i-1) = 1;
matrix(i,i) = A;
matrix(i,i+1) = 1;
rhs(i) = B;
end;
matrix(nunknowns, nunknowns-1) = 1; % last row
matrix(nunknowns, nunknowns) = A;
rhs(nunknowns) = B - Tb;
T = rhs/matrix; % solve for unknowns
Tfull(1) = Ta; % reconstruct full vector
Tfull(2:1 + nunknowns) = T(:);
Tfull(nunknowns + 2) = Tb;
position = 0:deltax:L;
write("position = ",position);
write("temperature distribution = ",Tfull);