function out = ex10ode(t,y,flag)

% Copyright 1999, The MathWorks, Inc.

global I P1 P3

if nargin < 3 | isempty(flag)
  G = 1;
  tau = 2.77e-4;
  xi = 0.8;
  A = 12.566;
  rho = 1000;
  V0 = 201.0619;
  R = 8314;
  C = 30354;
  mG = 8.397;
  Cv = 3.4153;
  g = 9.8;                % m/s^2  convert to m/h^2
  g = g*(60*60)^2;

  out = zeros(10,1);
  out(1) = (I*G - 2*xi*y(1)*tau - y(2))/tau^2;
  out(2) = y(1);
  out(3) = (y(7) - y(8))/A;
  out(4) = out(3)*A*y(6)*200/(mG*C);
  out(5) = y(5) - 0.03*exp(y(2)/0.28518);
  out(6) = y(7) - y(5)*Cv*sqrt(max(P1 - y(10),0));
  out(7) = y(8) - Cv*sqrt(max(y(10) - P3,0));
  out(8) = y(6) - mG*R*y(4)/(y(9)*1e3);
  out(9) = y(10) - y(6) - rho*g*y(3)/1e3; 
  out(10) = y(9) - V0 + A*y(3);
else
  switch(flag)           
  case 'mass'
    out = diag([1 1 1 1 0 0 0 0 0 0]);
  otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end       