function ex14
%                  DAE examples of index greater than 1
%
%  The pendulum problem arises most naturally as an index 3 problem
%  with an algebraic equation that says the length L of the pendulum
%  is constant.  It can be reduced by differentiation to an index 2
%  problem and by another to an index 1 problem.  The index 1 problem
%  is solved in ex1.  Here the higher index formulations are used to 
%  verify that dae will recognize a semi-explicit DAE of index greater
%  than 1 whilst finding consistent initial conditions.
%
%  The differential equations are
%
%        p' = u                                  
%        q' = v                                  
%       mu' = -p*lambda                          
%       mv' = -q*lambda - g   
%
%  The algebraic equation for index 3 is                   
%
%        0  = p^2 + q^2 - L^2 
%
%  The algebraic equation for index 2 is
%
%        0  = p*u + q*v   
%
%  The algebraic equation for index 1 is
%
%        0  = m*(u^2 + v^2) - g*q - lambda*L^2
%
%
%   Take  y(1) = p, y(2) = q, y(3) = u, y(4) = v, y(5) = lambda
%
%   Parameters    m = 1, L = 1, g = 1    and    0 <= t <= 10.

% Copyright 2004, The MathWorks, Inc.

lasterr('');

y0 = zeros(5,1);
y0(1) = 1;
y0(5) = 0.1;  %  Inconsistent initial condition.

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

% Problem parameter -- shared with the nested function
index = input('Select index 2 or 3:  ');

if index == 2 || index == 3
  try
    % This computation should not be succesful due to the index.
    [t y] = ode15s(@ex14ode,[0 10],y0,options);
  catch
    disp(lasterr);
    lasterr('');
  end
    
else
  fprintf('You have entered an index of %d, the solver has not been called.\n',index);
  fprintf('In this example, the solver will only be called with an index of 2 or 3.\n');
end

%------------------------------------------------------------------------------
% Nested function -- parameter 'index' shared with the outer function.

  function dy = ex14ode(t,y)

    yp = zeros(5,1);
    yp(1) =   y(3);
    yp(2) =   y(4);
    yp(3) = - y(1)*y(5);
    yp(4) = - y(2)*y(5) - 1;
    if index == 2
      yp(5) = y(1)*y(3) + y(2)*y(4);
    else
      yp(5) = y(1)^2 + y(2)^2 - 1;
    end
    dy = yp;
    
  end  % function ex14ode
  
%------------------------------------------------------------------------------

end  % ex14

