%  This function solves the fourth part of homework problem [H6.9]
%  ---------------------------------------------------------------
%
echo on
diary nsds_hw6_9d.dry
%
%  We use the standard Matlab stiff system solver
%
global sv st sw
sv = -1;
st = -1;
sw = +1;
tmax = 10;
[tvec,yvec] = ode15s('poiseuille3');
tvec = tvec';
yvec = yvec';
[n,m] = size(tvec);
out1 = poiseuille3(tmax,yvec(:,m));
err = norm(out1,'inf')
%
%  Plot the results
%
rmax = 1;
dr = rmax/50;
r = [0:dr:rmax]';
subplot(3,1,1)
plot(r,[yvec(50:-1:1,m);0])
grid on
title('ttt')
xlabel('xxx')
ylabel('yy1')
subplot(3,1,2)
plot(r,[0;yvec(150:-1:101,m)])
grid on
xlabel('xxx')
ylabel('yy2')
subplot(3,1,3)
plot(r,[yvec(100:-1:51,m);1])
grid on
xlabel('xxx')
ylabel('yy3')
print -deps nsds_hw6_9d.eps
%
diary off
echo off
%
return