function out = ex9ode(t,y,flag)

% Copyright 1999, The MathWorks, Inc.

if nargin < 3 | isempty(flag)
  B = [ 447.5 -452.5  -47.5  -52.5
       -452.5  447.5   52.5   47.5
        -47.5   52.5  447.5  452.5
        -52.5   47.5  452.5  447.5];

  B = B + 0.00025*[ 1  1  1  -1
                    1  1  1  -1
                    1  1  1  -1
                   -1 -1 -1   1];      
     
  out = zeros(8,1);
  r = (y(1) + y(2) + y(3) + y(4))/2;
  s = 0;
  for i = 1:4
    s = s + (r - y(i))^2;
  end  
  s = s/2;
  for i = 1:4
    out(i) = s - (r - y(i))^2;
  end
  out(1:4) = out(1:4) - B*y(1:4);
  out(5) = y(5) - y(1)*y(6);
  out(6) = 2*y(6) + y(6)^3 - y(1) + y(7) - 1 - exp(-t);
  out(7) = y(7) - y(8) + y(1)*y(6);
  out(8) = y(7) + y(8) + 5*y(1)*y(2);
else
  switch(flag)           
  case 'mass'
    out = diag([1 1 1 1 0 0 0 0]);
  otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end       