%
%  We use the standard Matlab stiff system solver
%
global sv st sw
sv = -1;
st = -1;
sw = +1;
tmax = 3.5;
[tvec,yvec] = ode15s('poiseuille2');

out3 = 

              AbsTol: []
                 BDF: []
              Events: []
         InitialStep: []
            Jacobian: []
           JConstant: []
            JPattern: []
                Mass: []
        MassConstant: []
        MassSingular: []
            MaxOrder: 2
             MaxStep: []
         NonNegative: []
         NormControl: 'on'
           OutputFcn: []
           OutputSel: []
              Refine: []
              RelTol: 1.0000e-005
               Stats: []
          Vectorized: []
    MStateDependence: []
           MvPattern: []
        InitialSlope: []

tvec = tvec';
yvec = yvec';
[n,m] = size(tvec);
out1 = poiseuille2(tmax,yvec(:,m));
err = norm(out1,'inf')

err =

    0.0500

%
%  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_9c.eps
%
diary off