clear
echo on
%In this example we show how a linear ODE
%may be solved using an implicit matrix method.
%The solution technique is based on making a
%finite difference approximation for the 
%derivatives, and then combining these to 
%convert the differential equation into a
%set of algebraic ones.  We do this with the
%following example out of heat transfer:
pause
%Consider a slab of width b in which a chemical
%reaction is going on.  The chemical reaction is
%exothermic, and is releasing energy at a rate
%S(x) which is a function of position.  This energy
%is diffusing out of the slab until it reaches
%the surface.  We assume both surfaces to be held
%at a constant temperature T0.  We wish to solve
%for the steady temperature distribution.
pause
%It is convenient to work in dimensionless variables
%(this makes any numerical solution to the problem
%more general), thus we define:
%  Tstar=(T-T0)/Tc
%  xstar=x/b
%  sstar=s/Sc
%where Sc is the characteristic magnitude of the source
%of energy, and finally:
%  Tc=Sc*b^2/k
%where k is the thermal conductivity.
pause
%The resulting differential equation is just:
%  Tstar'' = -sstar(xstar)
%with boundary conditions:
%  Tstar(0)=Tstar(1)=0
%We make use of the finite difference representation of
%this differential equation.  First, we specify the degree
%of discretization:
n=5;
dx=1/(n-1);
xstar=[0:dx:1]';
pause
%We now construct our matrix a:
a=-2*eye(n,n);
a=a+diag(ones(n-1,1),1);
a=a+diag(ones(n-1,1),-1);
%and finally divide by dx^2:
a=a/dx^2
pause
%The right hand side is easy.  We take the source
%to be just xstar^2:
sstar=-xstar.^2
pause
%To finish up, we use the boundary conditions at 
%zero and one.  This replaces the first and last
%rows of a and sstar:
a(1,1)=1;
a(1,2)=0;
a(n,n)=1;
a(n,n-1)=0;
sstar(1)=0;
sstar(n)=0;
pause
%The matrix a thus looks like:
a
%and the right hand side sstar is:
sstar
pause
%The solution is just:
Tstar=a\sstar
%or, plotting this up:
plot(xstar,Tstar,'o')
xlabel('position')
ylabel('Temperature')
pause
%In general, we want to do alot more discretization.
%We can simply increase n:
n=50;
%and we do the above all over again:
dx=1/(n-1);
xstar=[0:dx:1]';
a=-2*eye(n,n);
a=a+diag(ones(n-1,1),1);
a=a+diag(ones(n-1,1),-1);
a=a/dx^2;
sstar=-xstar.^2;
a(1,1)=1;a(1,2)=0;
a(n,n)=1;a(n,n-1)=0;
sstar(1)=0;sstar(n)=0;
Tstar=a\sstar;
%or, plotting this up:
hold on
plot(xstar,Tstar)
hold off
pause
%This gives us better resolution.  It is
%interesting to note the high degree of accuracy
%achieved with just a few points!  This is a
%consequence of the function being well behaved
%(no singularities in the higher order derivatives)
%and the O(dx^2) accuracy of the finite difference
%approximation.
pause
%This solution technique is very useful for
%linear ode's, and may be used for those where
%the coefficients are non-linear functions of x.
%Such equations are often difficult to solve
%analytically, and this provides a simple alternative
%solution approach.
echo off