daeic12: Need a better guess y0 for consistent initial conditions

I'm solving a system of DAEs over three intervals using ode15s as:

% Loop through each interval

for i = 1 : length(timePoints) - 1

    % Define the interval time span

    timeSpan = [timePoints(i) timePoints(i+1)] ;

    % Evaluate the initial slope and define the integration options
 
    dz = system.evaluateDAEs(timePoints(i),z0,u) ;
    options = odeset('Mass',M,'RelTol',1e-3,'AbsTol',1e-5,'MassSingular','yes','InitialSlope',dz,'MaxStep',1e-3) ;

    % Integrate over the current interval

    [tSol,zSol] = ode15s(@(t,z)system.evaluateDAEs(t,z,u),timeSpan,z0,options) ;

    % Update the solution

    tOutput = [tOutput
        tSol(2:end,:)] ;
    zOutput = [zOutput
        zSol(2:end,:)] ;

    % Extract the values of the state and algebraic variables at the end of
    % the current interval

    x = zOutput(end,1:numberOfStateVars) ;
    y = zOutput(end,numberOfStateVars+1:end) ;

    % Update the system to take into account potential disturbances at the
    % start of the new interval

    system = system.update(timePoints(i+1)) ; 

    % Solve for the algebraic variables after the disturbance

    [y,fval] = fsolve(@(y) system.evaluateAlgEqns(x,y,u,timePoints(i+1)),y) ;

    % Create a new initial vector

    z0 = [x y]' ;

end

Intervals are defined based on the occurrence of disturbances in the system. Disturbances are reflected as changes in the algebraic variable parameters. Since the state variables can't change immediately after a disturbance and algebraic variables can, I solve the algebraic equations using fsolve() each time after a disturbance to get new values of the algebraic variables. fsolve() solves the algebraic equations successfully. I use these values and the final values of the state variables from the previous interval to evaluate the initial slope used in ode15s. However, the solver returns an error:

Error using daeic12

Need a better guess y0 for consistent initial conditions.

Error in ode15s (line 298)

[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(ode,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Following the approach proposed here, I tried to get the initial slope as:

implicitDAE = @(yp) M*yp - system.evaluateDAEs(timePoints(i),z0,u) ;
dz0 = system.evaluateDAEs(timePoints(i),z0,u) ;
yp = fsolve(implicitDAE,dz0) ;

however, fsolve() just says that the equation is solved at the initial point. Can someone provide any insights how to solve this problem?

I should point out that this exact code works on a smaller system (lower number of state and algebraic equations).

[Lazar]

@Torsten, I saw that you have experience with this, so I would really appreciate your insight. Thanks in advance!

[Torsten]

If you don't supply executable code that reproduces the error (especially the functions in which algebraic and differential equations are defined), nobody will be able to help.