MATLAB: Solving a large system of DAEs using ode15s give error index greater than 1

MATLAB

When solving a large (N ~ 1000) system of Differential Algebraic Equations (DAEs) using the "ode15s" solver, I encountered an error "This DAE appears to be of index greater than 1". This error seems to be generated only when the value of a parameter "A" in the model equations is changed from 1e-10 to 1e-12. The system of equations is solved successfully when A = 1e-10, however when A=1e-12, the DAE error is generated.
Why does the index of equations become greater than 1 when the parameter value is changed?
Do you have any suggestions as to how to solve this issue ?
I am using the following options to solve my system of DAEs –
options = odeset('Mass',M,'MassSingular','yes','RelTol',1e-4,'AbsTol',1e6,'OutputFcn',@odetpbar);
[tout, yout] = ode15s(@f_pdes,tspan,y,options);

Best Answer

The “ode15s” solver requires the DAE to be of index 1.  The system is of index 1 when the matrix of partial derivatives of the algebraic system with respect to the algebraic variables is non-singular. The ode15s solver checks for this during the initialization phase and throws an error "This DAE appears to be of index greater than 1", as explained in this MATLAB Answers post- https://www.mathworks.com/matlabcentral/answers/102944-what-is-the-meaning-of-this-dae-appears-to-be-of-index-greater-than-1-using-ode-solvers-for-solvin  
In this case however, the DAE system is of index 1, since the ode15s is able to solve the system when the value of parameter A=1e-10. The ode15s solver may also generate this error when the "InitialSlope" used to initiate the solution is inaccurate. This issue can be addressed in a couple of ways:
1) Increasing the number of function evaluations using the "odeset" function:
When the "InitialSlope" is not provided as an option using the "odeset" function, the ode15s solver tries to compute an Initial Slope, which is required to advance the solution, using a Newton method. Increasing the number of "Maximum Function Evaluations" can improve the accuracy of the Initial Slope computed, especially for large systems. Instructions for setting the maximum number of function evaluations using the "odeset" function can be found in the following MATLAB Answers post - https://www.mathworks.com/matlabcentral/answers/94739-is-it-possible-to-specify-the-maximum-amount-of-time-or-number-of-function-evaluations-to-be-used-by
 
2) Providing a consistent "InitialSlope" to the ode15s solver:
An alternative approach is to provide the "InitialSlope" as an additional input parameter to the "ode15s" solver using the "odeset" function. The "fsolve" function in MATLAB, which solves a system of nonlinear equations, can be used to compute the initial slope. Moreover, the "fsolve" method offers greater control and flexibility in computing the InitialSlope using the "optimoptions" function.
For a system of the form:
M * (dy/dt) = systemDAE(t, y)
with the initial slope "yp" can be estimated by solving the non-linear equation
implicitDAE = @(yp) M*yp - systemDAE(0, y0)
where "y0" is a guess for the initial value of "y".
The initial slope can be computed using the "fsolve" function as:
yp = fsolve(implicitDAE, y0, options)
where "y0" is the initial guess provided for the initial slope, and the "optimoptions" function is used to set the options parameter.
Related Question