function  ex9
%  This is the DAS 1 test problem of I. Cameron, Solution of differential-
%  algebraic systems using diagonally implicit Runge-Kutta methods, IMA
%  J. Num. Anal., 3 (1983) 273-289.  He attributes it to C.W. Gear, 
%  Simultaneous numerical solution of differential-algebraic equations, 
%  IEEE Trans. Circuit Theory CT-18 (1971) 89-95.  The problem is to be 
%  solved on [0 1000].  Because the solution components tend quickly to 
%  limit values, they are plotted on a much shorter interval.

% Copyright 2002, The MathWorks, Inc.

%  Consistent initial conditions.
y0 = - ones(8,1);
y0(6) =  1;
y0(7) = -2;
y0(8) = -3;

M = diag([1 1 1 1 0 0 0 0]);
options = odeset('Mass',M);

[t,y] = ode15s(@ex9ode,[0 1000],y0,options);

i = find(t <= 6);
tlast = t(max(i));
figure(1)
plot(t(i),y(i,:))
axis([0 tlast -80 10])
title('Gear''s Stiff Test Problem')

%-------------------------------------------------------------------------%

function yp = ex9ode(t,y)

  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];      
     
  yp = 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
    yp(i) = s - (r - y(i))^2;
  end
  yp(1:4) = yp(1:4) - B*y(1:4);
  yp(5) = y(5) - y(1)*y(6);
  yp(6) = 2*y(6) + y(6)^3 - y(1) + y(7) - 1 - exp(-t);
  yp(7) = y(7) - y(8) + y(1)*y(6);
  yp(8) = y(7) + y(8) + 5*y(1)*y(2);
  