function out = ex13ode(t,y,flag)

% Copyright 1999, The MathWorks, Inc.

if nargin < 3 | isempty(flag)

  %  k = y(1), c = y(2), M = y(3), p = y(4), P = y(5), m = y(6)

  mu = 15 + 5*tanh(t - 10);
  mup = 5/cosh(t - 10)^2;

  f = zeros(6,1);
  f(1) =   (y(2) - y(1))/20;
  f(2) = - (y(4) - 99.1)/75;
  f(3) =    mu - y(6);
  f(4) =   (3.35 - 0.075*y(6) + 0.001*y(6)^2) - y(4)/y(5);
  f(5) =    20*y(5) - y(3);
  f(6) =    y(3)*f(3)/400 + mu*mup/(1.2*y(1))^2 ...
     - mu^2 * f(1)/(1.2^2 * y(1)^3);
  out = f;
else
  switch(flag)           
  case 'mass'
    M = diag([1 1 1 0 0 0]);
    M(2,4) = 1/15;
    out = M;
  otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end     