clear
reset(0)
echo on
%Example 31a
%
%Numerical stability problems often occur in the
%solution of partial differential equations.  Here
%we illustrate this with the example of the quenching
%of a finite slab using a very simple numerical method.
%
%Consider a slab of width 2b initially at a temperature T1,
%and then quench the temperature at the surface to a value
%T0.  The thermal diffusivity alpha is assumed to be a
%constant.  We render x dimensionless with respect to the
%half-width b, and time dimensionless with respect to the
%diffusion time b^2/alpha.  Thus we have the equation:
%
% dT/dt = d2T/dx2; T(t=0)=1, T(x=1)=0, dT/dx(x=0)=0
%
%The last is the symmetry condition at the centerline.
pause
%OK, now we solve this as a system of first order ODE's.
%We discretize x:
n=10;
dx=1/n;
x=[0:dx:1];

%and we also discretize T:
T=ones(size(x));
pause
%Now we solve this as a function of time.  We shall use
%the Euler method for simplicity:
t=0;
dt=0.005;
T(n+1)=0;
figure(1)
plot(x,T)
drawnow
for k=1:1/dt
 i=[2:n];  %We only apply the differential eqn to interior pts
 dT(i)=dt*(T(i-1)+T(i+1)-2*T(i))/dx^2;
 T(i)=T(i)+dT(i);
 T(n+1)=0;             %This is the zero Temperature BC at x=1
 T(1)=(4*T(2)-T(3))/3; %This is the zero derivative BC at x=0
 plot(x,T)
 axis([0 1 0 1])
 xlabel('x')
 ylabel('T')
 t=t+dt;
 drawnow
 echo off
end
echo on
pause
%The integration procedure behaved very well.  Now let's do
%it again with a -very- slightly larger time step:
dt=0.0052;
pause
t=0;
T=ones(size(x));
T(n+1)=0;
figure(1)
plot(x,T)
drawnow
for k=1:1/dt
 i=[2:n];
 dT(i)=dt*(T(i-1)+T(i+1)-2*T(i))/dx^2;
 T(i)=T(i)+dT(i);
 T(n+1)=0;
 T(1)=(4*T(2)-T(3))/3;
 plot(x,T)
 axis([0 1 0 1])
 xlabel('x')
 ylabel('T')
 t=t+dt;
 drawnow
 echo off
end
echo on
pause
%The dramatic difference between these two results illustrates
%the problem of numerical instability!  We require dt < 0.5 dx^2
%for the method to be stable.  Any explicit integration method
%would have the same stability limitations.  We can get around 
%this limitation by using implicit integration techniques.
echo off