% 1-D steady state heat transfer
% model d2T/dx2 = 0
% Dirichlet and Robin initial conditions 
% 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

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 equation is
A(1,1) = 1;
RHS(1) = BT0;
% for the Robin boundary conditions suppose we want the heat flux at the boundary to be proportional
% to the difference between T at boundary and the 'outside' temperature at the boundary
% , that is, dT/dx = -a*(T - Toutside)
% so a*T + dT/dx = a*Toutside
a = 1;
Toutside = 10;
A(11,11) = a + 1 / dx;
A(11,10) = -1 / dx;
RHS(11) = a*Toutside;  % the desired value of a*T + dT/dx so that dT/dx = -a*(T - Toutside)

% T at the interior points satisfied the DE d2T/dx2 = 0
% 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
% gives (T(i+1) - 2T(i) + T(i-1))/dx*dx = 0
for i=2:10
    A(i,i+1) = 1 / dx^2;
    A(i,i) = -2 / dx^2;
    A(i,i-1) = 1 / dx^2;
                         % RHS(i) = 0
end

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