Error for initial conditions for DAE Mass Matrix Solver

My mass matrix solver is being used to solve the states of piston engine. Due to exhaust and intake occuring during different crank angles I have to interrupt the ODE function to turn on and off the intake and exhaust, along with when the fuel reaction occurs. So I take the last values of the previous integration and input them as the initial conditions. This works for the most part but after several cycles, the ODE errors with:

Error using daeic3 (line 230)
Need a better guess y0 for consistent initial conditions.

Error in ode15s (line 313)
    [y,yp,f0,dfdy,nFE,nPD,Jfac,dMfac] = daeic3(odeFcn,odeArgs,tspan,htry,Mtype,Mt,Mfun,...

Error in IntegrationFunction (line 143)
[tSol,ySol,te,ye,ie] = ode15s(F, [ThetaStart, ThetaEnd], y0, opt);

Error in PistonEngineNonConstantCv_V4 (line 37)
    [tSol1,ySol1,S,origVars] = IntegrationFunction(ySol(a,1),ySol(a,2),ySol(a,3),...

There seems to be a trivial amount of cycles needed for this error to occur depending on the speed of the piston. The code can run for 6 cycles before erroring for 5 RPMs, run for 4 cycles for 10 RPMs, or will run for 7 cycles for 100 RPMs.

Here is my function for the DAE:

function [tSol,ySol,S,origVars] = IntegrationFunction(InitP,InitM,InitT,InitW,...
                       InitQ,InitV,InitS,ThetaStart,ThetaEnd,i)
%%This function takes the initial values from the Piston engine script
% and integrates the DAEs over the theta bounds provided

global RPM
syms p(theta) m(theta) T(theta) W(theta) Q(theta) V(theta) S(theta)...
        mdotin mdotex mdotleak Cp Tair Qdotrxn Qdotloss omega B L a

     %Air Properties
CylDia = .1;               %Cylinder diameter in m
R = 287.058/1000;          %Gas Constant for Air
BDCpres = 101324;          %Pressure of Cyclinder at BDC
PortDiameter = CylDia*(1/3);      %Intake and Exhaust port Diameter (m)
PortArea = pi*(PortDiameter^2)/4; %Intake and Exhaust Port Area (m)
C = 1;                     %Dischage Coefficient
Po = BDCpres;              %Upstream pressure            (N*m^2)
To = Tair;                 %Upstream temperature         (K)
gamma = 1.4;               %Specific Heat Ratio

if InitP <= Po
    Psi = real(sqrt(abs(gamma/(gamma-1)*((p(theta)/Po)^(2/gamma) -...
            (p(theta)/Po)^((gamma+1)/gamma)))));
else 
    Psi = real(sqrt(gamma/(gamma-1)*((Po/p(theta))^(2/gamma) -...
        (Po/p(theta))^((gamma+1)/gamma))));
end

mdotin = 0;
mdotex = 0;
mdotleak = 0;

if (ThetaStart == 0) && (ThetaEnd == pi)
    mdotin = C*PortArea*Psi*sqrt(2*Po*Po/(To*R*1000));
    mdotex = C*PortArea*Psi*sqrt(2*p(theta)*p(theta)/(R*1000*T(theta)));
    mdotleak = 0;
end

if (ThetaStart == pi*2-1/4*pi+pi*4*(i-1)) && (ThetaEnd == 2*pi+1/4*pi+pi*4*(i-1))
    Qdotrxn = 0;
else
    Qdotrxn = 0;
end

if  (ThetaStart == 3*pi) && (ThetaEnd == 4*pi)
    mdotin = C*PortArea*Psi*sqrt(2*Po*Po/(To*R*1000));
    mdotex = C*PortArea*Psi*sqrt(2*p(theta)*p(theta)/(R*1000*T(theta)));
    mdotleak = 0;
end 

%Constants of Cv Polynomial for Air
Alpha = 3.653;
Beta = -1.337*10^-3;
Gamma = 3.294*10^-6;
Delta = -1.913*10^-9;
Epsilon = 0.2763*10^-12;
CvT = R*(Alpha-1+Beta*T(theta)+ Gamma*T(theta)^2+Delta*T(theta)^3 +...
        Epsilon*T(theta)^4);
CpT = Alpha+Beta*T(theta)+ Gamma*T(theta)^2+Delta*T(theta)^3 +...
        Epsilon*T(theta)^4;
dCvT = 0;

    %Equation of State
eqn1 = diff(p(theta),1)/p(theta)-diff(m(theta), 1)/m(theta)...
        -diff(T(theta),1)/T(theta)+diff(V(theta), 1)/V(theta) == 0;
    %Conservation of Mass
eqn2 = diff(m(theta),1) == 1/omega*(mdotin - mdotex - mdotleak);
    %Conservation of Energy
eqn3 = CvT*T(theta)*diff(m(theta),1) +...
        m(theta)*(diff(T(theta),1)*CvT + T(theta)*dCvT) -...
        diff(Q(theta), 1)+diff(W(theta), 1) == 0;
    %Equation for Work
eqn4 = diff(W(theta),1)-p(theta)*diff(V(theta), 1) == 0;
    %Equation for Heat Exchange
eqn5 = diff(Q(theta),1) == 1/omega*(mdotin*CpT*Tair -...
        mdotex*CpT*T(theta) - mdotleak*CpT*T(theta) + Qdotrxn - Qdotloss);
    %Equation for Volume of Engine
eqn6 = diff(V(theta),1) == (B^2*pi*(a*sin(theta) +...
        (a^2*cos(theta)*sin(theta))/(- a^2*sin(theta)^2 + L^2)^(1/2)))/4;
    %Second Law of Thermodynamics
eqn7 = diff(S(theta),1) == diff(Q(theta),1)/T(theta);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7];

vars = [p(theta); m(theta); T(theta); W(theta); Q(theta); V(theta);...
        S(theta)];
origVars = length(vars);

[DAEs,DAEvars] = reduceDAEIndex(eqns,vars);

[M,f] = massMatrixForm(DAEs,DAEvars);

pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars);

M = odeFunction(M, DAEvars);
f = odeFunction(f, DAEvars, Tair ,...
                     Qdotloss, omega, B, L, a);
                
%Cyclinder Properties             
B = .1;          %Cylinder diameter in m
L = .15;         %Length of crank throw in m
a = .1250;          %Length of connecting rod in m
CompRatio = 10;
Vol_0 = pi*B*B*L/(4*(CompRatio-1));     %cylinder head volume       (m^3)
VolBDC = Vol_0+(pi*B*B*L/4);            %cylinder volume at BDC     (m^3)

%Air Properties
Tair = 300;                %Temperature of air in Kelvin
BDCpres = 101324;          %Pressure of Cyclinder at BDC

%Speed of Engine
rpm = RPM;             % Revolutions per a minute
RevTime = 60/rpm;       % Time for one revolution (in seconds)
omega = 2*pi/RevTime;   % Changing theta vs. time

%Heat exchange values
Qdotloss = 0;
QReleaseWidth = .1;
QReleaseTime = QReleaseWidth*RevTime;

F = @(theta, Y) f(theta, Y, Tair,...
              Qdotloss, omega, B, L, a);

% y0 = [round(InitP,5,'significant'); round(InitM,5,'significant');...
%         round(InitT,5,'significant'); round(InitW,5,'significant');...
%         round(InitQ,5,'significant'); round(InitV,5,'significant');...
%         round(InitS,5,'significant')];
          
y0 = [InitP; InitM; InitT; InitW; InitQ; InitV; InitS];
yp0 = zeros(7,1);


opt = odeset('Mass', M, 'InitialSlope', yp0,...
                'RelTol', 10.0^(-11), 'AbsTol', 10.0^(-11));
            
implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y);
% [y0, yp0] = decic(implicitDAE, 0.005, y0est, [], yp0est, [], opt);
            
opt = odeset(opt, 'InitialSlope', yp0, 'Events',...
            @PsiEventFunction);
            
[tSol,ySol,te,ye,ie] = ode15s(F, [ThetaStart, ThetaEnd], y0, opt);

x = length(ySol);
X = length(ySol)-1;
teTF = isempty(te);
timeevent1 = te*(180/pi)
te1 = [];
ySol(x,:) = real(ySol(x,:));

%%Post-first event if event happens
if teTF == 0
    syms p(theta) m(theta) T(theta) W(theta) Q(theta) V(theta) S(theta)...
        mdotin mdotex mdotleak Cp Tair Qdotrxn Qdotloss omega B L a

     %Air Properties
    CylDia = .1;               %Cylinder diameter in m
    R = 287.058/1000;          %Gas Constant for Air
    BDCpres = 101324;          %Pressure of Cyclinder at BDC
    PortDiameter = CylDia*(1/3);       %Intake and Exhaust port Diameter (m)
    PortArea = pi*(PortDiameter^2)/4; %Intake and Exhaust Port Area (m)
    C = 1;                     %Dischage Coefficient
    Po = BDCpres;              %Upstream pressure            (N*m^2)
    To = Tair;                 %Upstream temperature         (K)
    gamma = 1.4;               %Specific Heat Ratio
    
    
    if ySol(X,1) > Po
        ySol(x,1) = ySol(x,1)-ySol(x,1)/10000;
        Psi = real(sqrt(gamma/(gamma-1)*((p(theta)/Po)^(2/gamma) -...
            (p(theta)/Po)^((gamma+1)/gamma))));
    else 
        ySol(x,1) = ySol(x,1)+ySol(x,1)/10000;
        Psi = real(sqrt(gamma/(gamma-1)*((Po/p(theta))^(2/gamma) -...
            (Po/p(theta))^((gamma+1)/gamma))));
    end
    
    mdotin = 0;
    mdotex = 0;
    mdotleak = 0;

   if (ThetaStart == 0) && (ThetaEnd == pi)
        mdotin = C*PortArea*Psi*sqrt(2*Po*Po/(To*R*1000));
        mdotex = C*PortArea*Psi*sqrt(2*p(theta)*p(theta)/(R*1000*T(theta)));
        mdotleak = 0;
   end

   if (ThetaStart == pi*2-1/4*pi+pi*4*(i-1)) && (ThetaEnd == 2*pi+1/4*pi+pi*4*(i-1))
        Qdotrxn = 0;
   else
        Qdotrxn = 0;
   end

   if  (ThetaStart == 3*pi) && (ThetaEnd == 4*pi)
        mdotin = C*PortArea*Psi*sqrt(2*Po*Po/(To*R*1000));
        mdotex = C*PortArea*Psi*sqrt(2*p(theta)*p(theta)/(R*1000*T(theta)));
        mdotleak = 0;
   end 

    %Constants of Cv Polynomial for Air
    Alpha = 3.653;
    Beta = -1.337*10^-3;
    Gamma = 3.294*10^-6;
    Delta = -1.913*10^-9;
    Epsilon = 0.2763*10^-12;
    CvT = R*(Alpha-1+Beta*T(theta)+ Gamma*T(theta)^2+Delta*T(theta)^3 +...
            Epsilon*T(theta)^4);
    CpT = Alpha+Beta*T(theta)+ Gamma*T(theta)^2+Delta*T(theta)^3 +...
            Epsilon*T(theta)^4;
    dCvT = 0;

        %Equation of State
    eqn1 = diff(p(theta),1)/p(theta)-diff(m(theta), 1)/m(theta)...
        -diff(T(theta),1)/T(theta)+diff(V(theta), 1)/V(theta) == 0;
        %Conservation of Mass
    eqn2 = diff(m(theta),1) == 1/omega*(mdotin - mdotex - mdotleak);
        %Conservation of Energy
    eqn3 = CvT*T(theta)*diff(m(theta),1) +...
        m(theta)*(diff(T(theta),1)*CvT + T(theta)*dCvT) -...
        diff(Q(theta), 1)+diff(W(theta), 1) == 0;
        %Equation for Work
    eqn4 = diff(W(theta),1)-p(theta)*diff(V(theta), 1) == 0;
        %Equation for Heat Exchange
    eqn5 = diff(Q(theta),1) == 1/omega*(mdotin*CpT*Tair -...
        mdotex*CpT*T(theta) - mdotleak*CpT*T(theta) + Qdotrxn - Qdotloss);
        %Equation for Volume of Engine
    eqn6 = diff(V(theta),1) == (B^2*pi*(a*sin(theta) +...
        (a^2*cos(theta)*sin(theta))/(- a^2*sin(theta)^2 + L^2)^(1/2)))/4;
        %Second Law of Thermodynamics
    eqn7 = diff(S(theta),1) == diff(Q(theta),1)/T(theta);
    eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7];

    vars = [p(theta); m(theta); T(theta); W(theta); Q(theta); V(theta);...
            S(theta)];
    origVars = length(vars);

    [DAEs,DAEvars] = reduceDAEIndex(eqns,vars);

    [M,f] = massMatrixForm(DAEs,DAEvars);

    pDAEs = symvar(DAEs);
    pDAEvars = symvar(DAEvars);
    extraParams = setdiff(pDAEs,pDAEvars);

    M = odeFunction(M, DAEvars);
    f = odeFunction(f, DAEvars, Tair ,...
                     Qdotloss, omega, B, L, a);
                
    %Cyclinder Properties             
    B = .1;          %Cylinder diameter in m
    L = .15;         %Length of crank throw in m
    a = .1250;          %Length of connecting rod in m
    CompRatio = 10;
    Vol_0 = pi*B*B*L/(4*(CompRatio-1));     %cylinder head volume       (m^3)
    VolBDC = Vol_0+(pi*B*B*L/4);            %cylinder volume at BDC     (m^3)

    %Air Properties
    Tair = 300;                %Temperature of air in Kelvin
    BDCpres = 101324;          %Pressure of Cyclinder at BDC

    %Speed of Engine
    rpm = RPM;             % Revolutions per a minute
    RevTime = 60/rpm;       % Time for one revolution (in seconds)
    omega = 2*pi/RevTime;   % Changing theta vs. time

    %Heat exchange values
    Qdotloss = 0;
    QReleaseWidth = .1;
    QReleaseTime = QReleaseWidth*RevTime;

    F = @(theta, Y) f(theta, Y, Tair,...
              Qdotloss, omega, B, L, a);

    y01 = [ySol(x,1); ySol(x,2); ySol(x,3); ySol(x,4); ySol(x,5);...
            ySol(x,6); ySol(x,7)];
    yp0 = zeros(7,1);

    opt = odeset('Mass', M, 'InitialSlope', yp0,...
                'RelTol', 10.0^(-11), 'AbsTol', 10.0^(-11));
            
    implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y);
    % [y0, yp0] = decic(implicitDAE, 0.005, y0est, [], yp0est, [], opt);
            
    opt = odeset(opt, 'InitialSlope', yp0, 'Events',...
            @PsiEventFunction);
    
    [tSol1,ySol1,te1,ye1,ie1] = ode15s(F, [te, ThetaEnd], y01, opt);
    
    tSol = cat(1,tSol,tSol1);
    ySol = cat(1,ySol,ySol1);
    timeevent2 = te1*(180/pi)
end

%%Post-second event if event happens
teTF1 = isempty(te1);
if teTF1 == 0 
    x1 = length(ySol);
    X1 = length(ySol)-1;
    
        syms p(theta) m(theta) T(theta) W(theta) Q(theta) V(theta) S(theta)...
        mdotin mdotex mdotleak Cp Tair Qdotrxn Qdotloss omega B L a

     %Air Properties
    CylDia = .1;               %Cylinder diameter in m
    R = 287.058/1000;          %Gas Constant for Air
    BDCpres = 10134;          %Pressure of Cyclinder at BDC
    PortDiameter = CylDia*(1/3);       %Intake and Exhaust port Diameter (m)
    PortArea = pi*(PortDiameter^2)/4; %Intake and Exhaust Port Area (m)
    C = 1;                     %Dischage Coefficient
    Po = BDCpres;              %Upstream pressure            (N*m^2)
    To = Tair;                 %Upstream temperature         (K)
    gamma = 1.4;               %Specific Heat Ratio
    
    if ySol(X1,1) > Po
        ySol(x1,1) = ySol(x1,1)-ySol(x1,1)/10000;
        Psi = real(sqrt(gamma/(gamma-1)*((p(theta)/Po)^(2/gamma) -...
            (p(theta)/Po)^((gamma+1)/gamma))));
    else
        ySol(x1,1) = ySol(x1,1)+ySol(x1,1)/10000;
        Psi = real(sqrt(gamma/(gamma-1)*((Po/p(theta))^(2/gamma) -...
            (Po/p(theta))^((gamma+1)/gamma))));
    end
    
    mdotin = 0;
    mdotex = 0;
    mdotleak = 0;

    if (ThetaStart == 0) && (ThetaEnd == pi)
        mdotin = C*PortArea*Psi*sqrt(2*Po*Po/(To*R*1000));
        mdotex = C*PortArea*Psi*sqrt(2*p(theta)*p(theta)/(R*1000*T(theta)));
        mdotleak = 0;
    end

    if (ThetaStart == pi*2-1/4*pi+pi*4*(i-1)) && (ThetaEnd == 2*pi+1/4*pi+pi*4*(i-1))
        Qdotrxn = 0;
    else
        Qdotrxn = 0;
    end

    if  (ThetaStart == 3*pi) && (ThetaEnd == 4*pi)
        mdotin = C*PortArea*Psi*sqrt(2*Po*Po/(To*R*1000));
        mdotex = C*PortArea*Psi*sqrt(2*p(theta)*p(theta)/(R*1000*T(theta)));
        mdotleak = 0;
    end  

    %Constants of Cv Polynomial for Air
    Alpha = 3.653;
    Beta = -1.337*10^-3;
    Gamma = 3.294*10^-6;
    Delta = -1.913*10^-9;
    Epsilon = 0.2763*10^-12;
    CvT = R*(Alpha-1+Beta*T(theta)+ Gamma*T(theta)^2+Delta*T(theta)^3 +...
            Epsilon*T(theta)^4);
    CpT = Alpha+Beta*T(theta)+ Gamma*T(theta)^2+Delta*T(theta)^3 +...
            Epsilon*T(theta)^4;
    dCvT = 0;

        %Equation of State
    eqn1 = diff(p(theta),1)/p(theta)-diff(m(theta), 1)/m(theta)...
        -diff(T(theta),1)/T(theta)+diff(V(theta), 1)/V(theta) == 0;
        %Conservation of Mass
    eqn2 = diff(m(theta),1) == 1/omega*(mdotin - mdotex - mdotleak);
        %Conservation of Energy
    eqn3 = CvT*T(theta)*diff(m(theta),1) +...
        m(theta)*(diff(T(theta),1)*CvT + T(theta)*dCvT) -...
        diff(Q(theta), 1)+diff(W(theta), 1) == 0;
        %Equation for Work
    eqn4 = diff(W(theta),1)-p(theta)*diff(V(theta), 1) == 0;
        %Equation for Heat Exchange
    eqn5 = diff(Q(theta),1) == 1/omega*(mdotin*CpT*Tair -...
        mdotex*CpT*T(theta) - mdotleak*CpT*T(theta) + Qdotrxn - Qdotloss);
        %Equation for Volume of Engine
    eqn6 = diff(V(theta),1) == (B^2*pi*(a*sin(theta) +...
        (a^2*cos(theta)*sin(theta))/(- a^2*sin(theta)^2 + L^2)^(1/2)))/4;
        %Second Law of Thermodynamics
    eqn7 = diff(S(theta),1) == diff(Q(theta),1)/T(theta);
    eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7];

    vars = [p(theta); m(theta); T(theta); W(theta); Q(theta); V(theta);...
            S(theta)];
    origVars = length(vars);

    [DAEs,DAEvars] = reduceDAEIndex(eqns,vars);

    [M,f] = massMatrixForm(DAEs,DAEvars);

    pDAEs = symvar(DAEs);
    pDAEvars = symvar(DAEvars);
    extraParams = setdiff(pDAEs,pDAEvars);

    M = odeFunction(M, DAEvars);
    f = odeFunction(f, DAEvars, Tair ,...
                     Qdotloss, omega, B, L, a);
                
    %Cyclinder Properties             
    B = .1;          %Cylinder diameter in m
    L = .15;         %Length of crank throw in m
    a = .1250;          %Length of connecting rod in m
    CompRatio = 10;
    Vol_0 = pi*B*B*L/(4*(CompRatio-1));     %cylinder head volume       (m^3)
    VolBDC = Vol_0+(pi*B*B*L/4);            %cylinder volume at BDC     (m^3)

    %Air Properties
    Tair = 300;                %Temperature of air in Kelvin
    BDCpres = 101324;          %Pressure of Cyclinder at BDC

    %Speed of Engine
    rpm = RPM;             % Revolutions per a minute
    RevTime = 60/rpm;       % Time for one revolution (in seconds)
    omega = 2*pi/RevTime;   % Changing theta vs. time

    %Heat exchange values
    Qdotloss = 0;
    QReleaseWidth = .1;
    QReleaseTime = QReleaseWidth*RevTime;

    F = @(theta, Y) f(theta, Y, Tair,...
              Qdotloss, omega, B, L, a);

    y02 = [ySol(x1,1); ySol(x1,2); ySol(x1,3); ySol(x1,4); ySol(x1,5);...
            ySol(x1,6); ySol(x1,7)];
    yp0 = zeros(7,1);

    opt = odeset('Mass', M, 'InitialSlope', yp0,...
                'RelTol', 10.0^(-11), 'AbsTol', 10.0^(-11));
            
    implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y);
    % [y0, yp0] = decic(implicitDAE, 0.005, y0est, [], yp0est, [], opt);
            
    opt = odeset(opt, 'InitialSlope', yp0, 'Events',...
            @PsiEventFunction);
        
    [tSol2,ySol2,te2,ye2,ie2] = ode15s(F, [te1, ThetaEnd], y02, opt);
    
    tSol = cat(1,tSol,tSol2);
    ySol = cat(1,ySol,ySol2);
    timeevent3 = te2*(180/pi)
end     

%%Organize ySol and tSol for plotting
S = cell(origVars, 1); %For labels of plots

for k = 1:origVars
  S{k} = char(DAEvars(k));
end

tSol = tSol.*(180/pi);

end

The reason for the multiple interations of the mass matrix solver is due to the coefficient Psi, which changes when pressure inside of the cylinder equals atmospheric pressure, or else Psi would become imaginary.

Here is the event function that interrupts the ODE solver:

function [val, isterminal, dir] = PsiPositiveEventFunction(theta,p)
%%This is an event function for ode15s to interrupt the integration
%when Po = p(theta) so that Psi does not become imaginary
Po = 101324;  %Atmospheric Pressure (N/m^2)


val = p(1) - Po;    %The value that we want to be zero for the event
isterminal = 1;     %This halts the integration if event is met
dir = 0;            %p(theta) approaches from both sides

And here is the code that runs the ODE function:

clear all;
global RPM

%Universal Gas Constant
R = 287.058/1000;
             
%Cyclinder Properties             
B = .1;          %Cylinder diameter in m
L = .15;         %Length of crank throw in m
a = .125;        %Length of connecting rod in m
V0 = .001;       %Clearance Height m
CompRatio = 10;
Vol_0 = pi*B*B*L/(4*(CompRatio-1));     %cylinder head volume       (m^3)
VolBDC = Vol_0+(pi*B*B*L/4);            %cylinder volume at BDC     (m^3)

%Air Properties
BDCpres = 101323;      %Atmospheric pressure (N/m^2)
Tair = 300;      %Temperature of air in Kelvin
initm = (Vol_0*BDCpres)/(Tair*R);

%Fuel Properties
fuel_AFR = 14.6;                % Diesel AFR
comb_eff = .95;                 % Assumed combustion efficiency
fuel_specific_energy = 48e6;    % Diesel specific energy, J/kg
fuel_energy = comb_eff*fuel_specific_energy; % Available specific energy

Cycles = 10;
RPM = 2500;

ySol = [BDCpres initm Tair 0 0 Vol_0 0];
tSol = zeros(1);

for i = 1:Cycles
    %%For 0<theta<180 with open intake, 1st Cycle
    a = length(ySol(:,1));
    
    [tSol1,ySol1,S,origVars] = IntegrationFunction(ySol(a,1),ySol(a,2),ySol(a,3),...
        ySol(a,4),ySol(a,5),ySol(a,6),ySol(a,7),0,pi,i);

    %%For 180<theta<315 with closed intake and exhaust
    a = length(ySol1); 
    
    [tSol2,ySol2,S,origVars] = IntegrationFunction(ySol1(a,1),ySol1(a,2),ySol1(a,3),...
        ySol1(a,4),ySol1(a,5),ySol1(a,6),ySol1(a,7),pi,pi*2-1/4*pi,i);
    
    %%For 315<theta<405 Ignition with closed intake and exhaust
    a = length(ySol2);
    
    [tSol3,ySol3,S,origVars] = IntegrationFunction(ySol2(a,1),ySol2(a,2),ySol2(a,3),...
        ySol2(a,4),ySol2(a,5),ySol2(a,6),ySol2(a,7),pi*2-1/4*pi,pi*2+1/4*pi,i);

    %%For 405<theta<540 with closed intake and exhaust
    a = length(ySol3); 
    
    [tSol4,ySol4,S,origVars] = IntegrationFunction(ySol3(a,1),ySol3(a,2),ySol3(a,3),...
        ySol3(a,4),ySol3(a,5),ySol3(a,6),ySol3(a,7),pi*2+1/4*pi,pi*3,i);

    %%For 540<theta<720 with open exhaust
    a = length(ySol4); 
    
    [tSol5,ySol5,S,origVars] = IntegrationFunction(ySol4(a,1),ySol4(a,2),ySol4(a,3),...
        ySol4(a,4),ySol4(a,5),ySol4(a,6),ySol4(a,7),pi*3,pi*4,i);
    
    
    tSol = cat(1,tSol,tSol1+720*(i-1),tSol2+720*(i-1),tSol3+720*(i-1),tSol4+720*(i-1),tSol5+720*(i-1));
    ySol = cat(1,ySol,ySol1,ySol2,ySol3,ySol4,ySol5);
end

PlotFunction(S,tSol,ySol,origVars); %Plot function to plot the 7 variables dependent on theta

I have no idea what is causing this error due to these values coming from a previous integration so the conditions should still apply.