function out = ex6ode(t,y,flag)

% Copyright 1999, The MathWorks, Inc.

global y0

if nargin < 3 | isempty(flag)
  out = zeros(3,1);
  out(1) = - y(2)^2;

  R = 45.72;
  rho = 0.814;
  mu = 0.098;
  b = 0.345;
  Qco = 1.7153e6;
  Qa = 3.027e5;

  temp1 = pi*sqrt(0.5*R/rho);
  temp2 = sqrt(0.5*R*rho);
  temp3= log(temp2 * (y(3)/mu) * y(2) - 5);

  out(2) =   temp1*(R - y(3))^2 * y(2) * (2.5*temp3 + 10.5) ...
           - b*Qco - (y0(1)/y(1)) * Qco * (1 - b);

  out(3) = 2*temp1*y(2)*((2.5*R*y(3) - 1.25*y(3)^2)* ...
           temp3 + 3*R*y(3) - 2.125*y(3)^2 - 13.6*R*mu/(temp2*y(2))) - Qa;
else
  switch(flag)           
  case 'mass'
    out = diag([1 0 0]);
  otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end     