function [out1,out2,out3] = poiseuille2(tt,x,flag)
%
%  State-space model of second part of [H6.9]
%
global sv st sw
if nargin < 3 | isempty(flag)
  gama = 10;
  smax = 1;
  ds = smax/50;
  s = [0:ds:smax]';
  r = ones(size(s)) - s;
  v(2:51) = x(1:50);
  t(2:51) = x(51:100);
  w(1:50) = x(101:150);
  v(1) = 0;
  t(1) = 1;
  w(51) = 0;
  h1 = t(51) + 1;
  h2 = h1*h1;
  h3 = h2*h1;
  cv = sqrt(2*gama)/h2;
  cw = gama/h3;
  dw(1) = sw*((-w(3) + 4*w(2) - 3*w(1))/(2*ds) - cw*r(1)*r(1)*r(1)*t(1)*t(1));
  for i=2:50,
    av = r(i)*t(i)*t(i);
    aw = av*r(i)*r(i);
    dv(i) = sv*((v(i+1) - v(i-1))/(2*ds) - cv*av);
    dt(i) = st*((t(i+1) - t(i-1))/(2*ds) + t(i)*w(i)/r(i));
    dw(i) = sw*((w(i+1) - w(i-1))/(2*ds) - cw*aw);
  end
  dv(51) = sv*((3*v(51) - 4*v(50) + v(49))/(2*ds));
  dt(51) = st*((3*t(51) - 4*t(50) + t(49))/(2*ds));
  out1 = [dv(2:51)' ; dt(2:51)' ; dw(1:50)'];
 else
  switch(flag)
   case('init')
    tmax = 3.5;
    smax = 1;
    ds = smax/50;
    out1 = [0,tmax];
    out2 = [ds:ds:2,-1:ds:-ds]';
    out3 = odeset('RelTol',1.0e-5,'MaxOrder',2,'NormControl','on')
   case('jacobian')
    gama = 10;
    smax = 1;
    ds = smax/50;
    s = [0:ds:smax]';
    r = ones(size(s)) - s;
    v(2:51) = x(1:50);
    t(2:51) = x(51:100);
    w(1:50) = x(101:150);
    v(1) = 0;
    t(1) = 1;
    w(51) = 0;
    h1 = t(51) + 1;
    h2 = h1*h1;
    h3 = h2*h1;
    cv = sqrt(2*gama)/h2;
    cw = gama/h3;
    d = 1/(2*ds);
    out1 = zeros(150,150);
    out1(1,2) = sv*d;
    out1(1,51) = -sv*2*cv*r(2)*t(2);
    for i=2:49,
      out1(i,i-1) = -sv*d;
      out1(i,i+1) = sv*d;
      out1(i,i+50) = -sv*2*cv*r(i+1)*t(i+1);
    end
    out1(50,48) = sv*d;
    out1(50,49) = -sv*4*d;
    out1(50,50) = sv*3*d;
    out1(51,51) = st*w(2)/r(2);
    out1(51,52) = st*d;
    out1(51,102) = st*t(2)/r(2);
    for i=2:49,
      ii = i + 50;
      out1(ii,ii-1) = -st*d;
      out1(ii,ii) = st*w(i+1)/r(i+1);
      out1(ii,ii+1) = st*d;
      out1(ii,ii+51) = st*t(i+1)/r(i+1);
    end
    out1(100,98) = st*d;
    out1(100,99) = -st*4*d;
    out1(100,100) = st*3*d;
    out1(101,50) = -sw*2*cw*r(1)*r(1)*r(1)*t(1);
    out1(101,101) = -sw*3*d;
    out1(101,102) = sw*4*d;
    out1(101,103) = -sw*d;
    for i=2:49,
      ii = i + 100;
      out1(ii,ii-51) = -sw*2*cw*r(i)*r(i)*r(i)*t(i);
      out1(ii,ii-1) = -sw*d;
      out1(ii,ii+1) = sw*d;
    end
    out1(150,99) = -sw*2*cw*r(50)*r(50)*r(50)*t(50);
    out1(150,149) = -sw*d;
   otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end
%
return