function [out1,out2,out3] = poiseuille4(tt,x,flag)
%
%  State-space model of fourth part of [H6.9]
%
%
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) = (-w(3) + 4*w(2) - 3*w(1))/(2*ds) - cw*r(1)*r(1)*r(1)*t(1)*t(1);
  dv(2) = cv*r(2)*t(2)*t(2) - (v(3) - v(1))/(2*ds);
  dt(2) = -(t(3) - t(1))/(2*ds) - t(2)*w(2)/r(2);
  dw(2) = (-w(4) + 4*w(3) - 3*w(2))/(2*ds) - cw*r(2)*r(2)*r(2)*t(2)*t(2);
  for i=3:49,
    av = r(i)*t(i)*t(i);
    aw = av*r(i)*r(i);
    dv(i) = cv*av - (3*v(i) - 4*v(i-1) + v(i-2))/(2*ds);
    dt(i) = -(3*t(i) - 4*t(i-1) + t(i-2))/(2*ds) - t(i)*w(i)/r(i);
    dw(i) = (-w(i+2) + 4*w(i+1) - 3*w(i))/(2*ds) - cw*aw;
  end
  dv(50) = cv*r(50)*t(50)*t(50) - (3*v(50) - 4*v(49) + v(48))/(2*ds);
  dt(50) = -(3*t(50) - 4*t(49) + t(48))/(2*ds) - t(i)*w(i)/r(i);
  dw(50) = (w(51) - w(49))/(2*ds) - cw*r(50)*r(50)*r(50)*t(50)*t(50);
  dv(51) = -(3*v(51) - 4*v(50) + v(49))/(2*ds);
  dt(51) = -(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 = 10;
    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) = -d;
    out1(1,51) = 2*cv*r(2)*t(2);
    out1(2,1) = 4*d;
    out1(2,2) = -3*d;
    out1(2,52) = 2*cv*r(3)*t(3);
    for i=3:49,
      out1(i,i-2) = -d;
      out1(i,i-1) = 4*d;
      out1(i,i+1) = -3*d;
      out1(i,i+50) = 2*cv*r(i+1)*t(i+1);
    end
    out1(50,48) = -d;
    out1(50,49) = 4*d;
    out1(50,50) = -3*d;
    out1(51,51) = -w(2)/r(2);
    out1(51,52) = -d;
    out1(51,102) = -t(2)/r(2);
    out1(52,51) = 4*d;
    out1(52,52) = -3*d - w(3)/r(3);
    out1(52,103) = -t(3)/r(3);
    for i=3:49,
      ii = i + 50;
      out1(ii,ii-2) = -d;
      out1(ii,ii-1) = 4*d;
      out1(ii,ii) = -3*d - w(i+1)/r(i+1);
      out1(ii,ii+51) = -t(i+1)/r(i+1);
    end
    out1(100,98) = -d;
    out1(100,99) = 4*d;
    out1(100,100) = -3*d;
    out1(101,50) = -2*cw*r(1)*r(1)*r(1)*t(1);
    out1(101,101) = -3*d;
    out1(101,102) = 4*d;
    out1(101,103) = -d;
    for i=2:48,
      ii = i + 100;
      out1(ii,ii-51) = -2*cw*r(i)*r(i)*r(i)*t(i);
      out1(ii,ii) = -3*d;
      out1(ii,ii+1) = 4*d;
      out1(ii,ii+2) = -d;
    end
    out1(149,98) = -2*cw*r(49)*r(49)*r(49)*t(49);
    out1(149,149) = -3*d;
    out1(149,150) = 4*d;
    out1(150,99) = -2*cw*r(50)*r(50)*r(50)*t(50);
    out1(150,149) = -d;
   otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end
%
return