clear
reset(0)
echo on
%Solution to Problem 31
%
%Here we adapt example 31a to solve for startup flow in a
%channel of width 2b.  The problem is virtually identical
%to the relaxation of an initial temperature distribution
%in the same geometry, except that now we have to add in a
%source term.  We solve this by discretizing the domain in
%the y direction, and then develop a system of first order
%equations for velocity as a function of t.  We shall use
%the symmetry condition at the centerline and impose a zero
%derivative boundary condition at y=0.

pause
%OK, now we solve this as a system of first order ODE's.
%We discretize y:
n=10;
dy=1/n;
y=[0:dy:1];

%and we also discretize u:
u=zeros(n+1,1); %This is the initial velocity distribution
pause
%Now we solve this as a function of t.  We shall use
%the Euler method for simplicity:
%
%We construct the tri-diagonal matrix governing the
%discretized equation:
a=-2*eye(n+1);
a=a+diag(ones(1,n),-1);
a=a+diag(ones(1,n),+1);
a=a/dy/dy;
%And we modify the first and last rows to account for the
%boundary conditions:
a(1,:)=zeros(1,n+1);
a(n+1,:)=zeros(1,n+1);
pause
%We can now integrate forward in t:
t=0;
dt=0.4/n^2;  %We pick a dt for stable integration
figure(1)
plot(y,u)
drawnow
%We also want to keep all the t's and u's for later plotting:
tall=zeros(1/dt+1,1);
uall=zeros(1/dt+1,length(u));
tall(1)=t;
uall(1,:)=u';  %We have to transpose u to make it a row vector.
for k=1:1/dt
 dudt=a*u+ones(size(u));  %This is the time derivative with source term.
 u=u+dt*dudt;  %We update the velocity at the interior points.
 u(n+1)=0;             %This is the zero velocity BC at y=1
 u(1)=(4*u(2)-u(3))/3; %This is the zero derivative BC at y=0
 plot(y,u)
 axis([0 1 0 .5])
 xlabel('y')
 ylabel('u')
 t=t+dt;
 drawnow
 tall(k+1)=t;
 uall(k+1,:)=u';
 echo off
end
echo on
pause

%Now we do some plotting of the velocity field at appropriate
%values of t.  Because the discretization in t is so small, we
%want to use only some of the values:
i=[1:1/10/dt:1/dt+1];
plot(tall(i),uall(i,:)') %Note that we have to use uall transpose
xlabel('y')
ylabel('u')
title('Velocity Profile for Startup Flow in a Channel')
%We also want to produce a legend:
labels=['t = ',num2str(tall(1))];
for k=2:length(i)
 labels=str2mat(labels,['t = ',num2str(tall(i(k)))]);
 echo off
end
echo on
legend(labels)
%And that's it!