%
% This solves a steady, 1-D heat conduction problem with a source term.
% The governing equations are:
%   \lambda \frac{d^2}{dT^2} = s
% with 
%   s = \exp\left( -\frac{\left(z-\frac{L}{2}\right)^2}{\gamma} \right)
%
% Author: James C. Sutherland
%

clear; clc; close all;

%---------------------------------------------------------------
% problem setup:
n      = 200;  % number of cells
L      = 1;    % domain length
lambda = 1e-5; % thermal diffusivity
gamma  = 0.1;% controls width of the source term
%---------------------------------------------------------------


%---------------------------------------------------------------
% set cell spacing and the cell centroids
% note that the cell centroids are offset from the domain boundary by dz/2.
dz = L/n;
z = (dz/2 : dz : L-dz/2)';  % column vector of mesh positions.
%---------------------------------------------------------------


%---------------------------------------------------------------
% populate the matrix.  We will form the main, lower, and upper
% diagonals and use these to build the matrix.
mainDiagonal = -2*ones(n,1);
lowerDiagonal = ones(n-1,1);
upperDiagonal = ones(n-1,1);

% adjust at boundaries
mainDiagonal(1) = -3;
mainDiagonal(n) = -1;

A = diag(mainDiagonal,0)   ...
  + diag(lowerDiagonal,-1) ...
  + diag(upperDiagonal,1);
%---------------------------------------------------------------


%---------------------------------------------------------------
% set the rhs and solve for temperature
s = exp( -(z-L/2).^2 / gamma );
rhs = -dz^2/lambda * s;
T = A\rhs;
%---------------------------------------------------------------


plot(z,T,'o-');
xlabel('Position');
ylabel('Temperature');