finitedf.m

%
% Finite Difference method example of solving boundary value ODE
%
%     d^2 T   
%     -----  + c (Tinf - T) = 0
%      dx^2
%
% which models the temperature distribution in a rod

% define constants

Ta   = 40;      % left side temperature
Tb   = 200;     % right side temperature
Tinf = 20;      % ambient temperature
c    = 0.01;    % constant

L      = 10;    % length of bar

ndivs  = input('Input number of divisions: ');
nunknowns = ndivs - 1;
deltax = L/ndivs;

A = -(2 + deltax^2*c);
B = -deltax^2*c*Tinf;

matrix = zeros(nunknowns);  % initialize

matrix(1,1) = A;                          % first row
matrix(1,2) = 1;
rhs(1)      = B - Ta;

for i = 2:nunknowns - 1                   % intermediate rows
  matrix(i,i-1) = 1;
  matrix(i,i)   = A;
  matrix(i,i+1) = 1;
  rhs(i)        = B;
end;

matrix(nunknowns, nunknowns-1) = 1;       % last row
matrix(nunknowns, nunknowns)   = A;
rhs(nunknowns)                 = B - Tb;

T = matrix\rhs';                          % solve for unknowns

Tfull(1)               = Ta;              % reconstruct full vector
Tfull(2:1 + nunknowns) = T(:);
Tfull(nunknowns + 2)   = Tb;

position = 0:deltax:L;
plot(position,Tfull);

% Tfull' = 
%   40.0000
%   65.9698
%   93.7785
%  124.5382
%  159.4795
%  200.0000