% 1-D steady state heat transfer
% model d2T/dx2 = 0
% Dirichlet boundary conditions - with added heat sources and heat disappation to outside air
% heat transfer through an insulated wire of lenght 1, with one end held at T0 degrees 
% and the other at T1 degrees
L = 7;
T = zeros(L,1);    % the estimated temperature values
RHS = zeros(L,1);  % right hand side of the point equations
A = zeros(L,L);   % the array of coefficients of the linear equations for T1, T2 ..   

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

dx = .5; % 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:3;

% 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(L,L) = 1;
RHS(L) = BT1;   

% for the interior points the T solves the PDE d2T/dx2 + a*(T - Tamb) = 0
% substituting the central difference approximation for the 2nd derivative of T, gives
% (T(i+1) - 2T(i) + T(i-1))/dx*dx - a*(T(i) - Tamb) = 0 that the equation for the node is
%  
a = 1;
Tamb = 100;
for i=2:L-1
    A(i,i+1) = 1 / dx^2;
    A(i,i) = -2 / dx^2 - a;
    A(i,i-1) = 1 / dx^2;
    RHS(i) = -a*Tamb;
end


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