% 1-D steady state heat transfer
% model d2T/dx2 = 0
% Dirichlet boundary conditions - with added heat sources
% heat transfer through an insulated wire of lenght 1, with one end held at T0 degrees 
% and the other at T1 degrees

T = zeros(11,1);    % the estimated temperature values
RHS = zeros(11,1);  % right hand side of the point equations
A = zeros(11,11);   % the array of coefficients of the linear equations for T1, T2 ..   

BT0 = 0; % boundary temperature at x = 0
BT1 = 0; % boundary temperature at x = 1

dx = .1; % the grid points are 0, .1, .2 .3 .4 .5 .6 .7 .8 .9 1
         % note that there are 11 points in the grid numbered 1 2 3 .... 11
xi=0:dx:1;

% now we create an array A of the coefficients of the equations defining T(i)
% for the essential Dirichlet boundary conditions the equations are
A(1,1) = 1;     
RHS(1) = BT0;
A(11,11) = 1;
RHS(11) = BT1;   

% substituting the central difference approximation for the 2nd derivative of T, that is
% d2T/dx2 = (T(i+1) - 2T(i) + T(i-1))/dx*dx
% this equation applies to the interior points of the grid
for i=2:10
    A(i,i+1) = 1 / dx^2;
    A(i,i) = -2 / dx^2;
    A(i,i-1) = 1 / dx^2;
end

% add heat sources at the 6th grid point (.5)
RHS(6) = -200;

% let MATLAB compute the solution
T = inv(A)*RHS;
plot(xi,T)
xlabel('x')
ylabel('T')