%-----------------------------------------------------------------------
      function [r, T] = T (r_in, r_out, T_in, T_out, n)
%-----------------------------------------------------------------------
% The following differential equation describes the steady-state temperature
% profile across the pipe wall set up by heat conduction.
%
%   dT/dr + r (d^2T/dr^2) = 0       BC: T (r_in)  = T_in
%                                       T (r_out) = T_out
%
% The finite difference representation is:
%
%   (1/2+r/dr)*T(i+1) - 2*r/dr*T(i) + (-1/2+r/dr)*T(i-1) = 0   i=2,3,...,10
%
% where r = r_in+dr*(i-1)       i=1, 2,..., n+1
%
% An example with n=10
% i=1:  T(1) = T_in        ... specified
% i=2:  (-0.5+0.055/0.005)*T(1) -2*0.055/0.005*T(2) + (0.5+0.055/0.005)*T(3) = 0
% i=3:  (-0.5+0.060/0.005)*T(2) -2*0.060/0.005*T(3) + (0.5+0.060/0.005)*T(4) = 0
% i=4:  (-0.5+0.065/0.005)*T(3) -2*0.065/0.005*T(4) + (0.5+0.065/0.005)*T(5) = 0
%  :                                                                         :
% i=10: (-0.5+0.095/0.005)*T(9) -2*0.095/0.005*T(10)+ (0.5+0.095/0.005)*T(11)= 0
% i=n:  T(n) = T_out       ... specified
%
% Variables used ...
%   r_in  = radius of the inner wall in meter (input)
%   r_out = radius of the outer wall in meter (input)
%   T_in  = temperature at the inner wall in Celsius (input)
%   T_out = temperature at the outer wall in Celsius (input)
%   n     = number of divisions or intervals (input)
%   r(i)  = radius at point i (output)
%   T(i)  = temperature at point i (output)
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Declare variables ----------------------------------------------------
      r = zeros(n+1,1);
%     T = zeros(n+1,1);
      A = zeros(n+1,n+1);
      b = zeros(n+1,1);

% Assign constants in the linear equation ------------------------------
      dr = (r_out-r_in)/n;
      for i=2: n
        r(i) = r_in + dr*(i-1);
        A(i,i-1) = -1./2. + r(i)/dr;
        A(i,i)   = -2.*r(i)/dr;
        A(i,i+1) =  1./2. + r(i)/dr;
      end

% B.C. -----------------------------------------------------------------
      r(1)       = r_in;
      A(1,1)     = 1.;
      b(1)       = T_in;
      r(n+1)     = r_out;
      A(n+1,n+1) = 1.;
      b(n+1)     = T_out;

% Solve the linear set of equations with predivision -------------------
      T = A\b;
% Equally valid ways of finding solution to a linear system of equations.
%     T = inv(A)*b;
%     T = A^(A)*b;