%  This is the pendulum problem posed as an index 1 problem.  Formulated as 
%  an index 2 problem, it is used as an example in E. Hairer, Ch. Lubich,
%  and M. Roche, The Numerical Solution of Differential-Algebraic Systems 
%  by Runge-Kutta Methods, Lecture Notes in Math. #1409, Springer, Berlin,
%  1989. The values of the parameters here are taken from their FORTRAN 
%  program.  The problem is solved in forms of index 3, 2, and 1 in 
%  K. Brenan, S. Campbell, and L. Petzold, Numerical Solution of Initial-
%  Value Problems in Differential-Algebraic Equations, SIAM, Philadelphia,
%  1996.
%
%  (p,q) are the Cartesian coordinates of the bob of the pendulum.  (u,v)
%  are the velocities.  L is the length of the pendulum and lambda is the
%  tension in the rod.  g is the acceleration due to gravity.
%
%        p' = u                                  
%        q' = v                                  
%       mu' = -p*lambda                          
%       mv' = -q*lambda - g                      
%        0  = m*(u^2 + v^2) - g*q - l^2 * lambda   algebraic equation
%
%  Here  y(1) = p, y(2) = q, y(3) = u, y(4) = v, y(5) = lambda
%
%  The parameters are    m = 1, l = 1, g = 1    and    0 <= t <= 10.

% Copyright 1999, The MathWorks, Inc.

%  Hairer & Wanner's RADAU5 requires consistent initial conditions,
%  which are obvious in this case.
y0 = zeros(5,1);
y0(1) = 1;
options = odeset('Mass','M','MassSingular','yes');

[t,y] = ode15s('ex3ode',[0 10],y0,options);

plot(t,y(:,1:4))
title('Pendulum Problem as Index 1 DAE')

