Integration tolerances using ode23s

Question

I have to solve a differential function using an ode solver (preferably ode23S) in order to plot the change of variables. Using fsolve earlier in the code to obtain initial (steady state) variables and then implementing a step-change in order to see the change of these variables. When trying to solve the ODE I get the following error: Warning: Failure at t=1.956374e+03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.950436e-12) at time t.

clear
Parameters
c3 = 1.5; % mol^(0.5).s^-1
TF = 115; % °C (Given as constant)
T3 = 115; % °C (Given as constant)
cp = 294.44; % J.mol^-1.K^-1
gammav = 34.37 * 10^3; % J/mol
n0 = 900; % mol
c2 = 2; % mol^0.5/s
c1 = 1.2; % mol^0.5/s
alpha = 50; % mol/s.K
T2sat = 116.5; % °C
T1sat = 117.6; % °C
T0sat = 120.9; % °C
LF = 20; % mol/s
Q_ss = 1 * 10^6; % W

Differential equations
% X(1) = T0
% X(2) = T1
% X(3) = T2
% X(4) = n1
% X(5) = n2
% X(6) = n3
dX_dt = @(t, X) [((((c1.*sqrt(X(4))).*X(2))-((c1.*sqrt(X(4)) - alpha.*(X(1)-T0sat)).*X(1))-((alpha.*(X(1)-T0sat)).*(X(1)+(gammav/cp)))+(Q_ss/cp))/900); ...
                 (c2.*sqrt(X(5)))-(c1.*sqrt(X(4)))-(alpha.*(X(2)-T1sat)); ...
                 (((c2.*sqrt(X(5)).*(X(3)-X(2)))-((alpha.*(X(2)-T1sat)).*(gammav/cp))))./(((c1*sqrt(X(4)))/c1)^2); ...
                 LF+((c3.*sqrt(X(6)))./2)-(c2.*sqrt(X(5)))+(alpha.*(X(1)-T0sat))+(alpha.*(X(2)-T1sat))-(alpha.*(X(3)-T2sat)); ...
                 ((LF.*(TF-X(3)))+(((c3.*sqrt(X(6)))./2)).*(T3-X(3))+((alpha.*(X(1)-T0sat)).*(X(1)-X(3)+(gammav/cp)))+((alpha.*(X(2)-T1sat)).*(X(2)-X(3)+(gammav/cp)))-((alpha.*(X(3)-T2sat)).*(gammav/cp)))./((((c2*sqrt(X(5)))/c2)^2)); ...
                 (alpha.*(X(3)-T2sat))-2*((c3.*sqrt(X(6)))./2)];

 Steady state values
X_initial = fsolve(@(X) dX_dt(0, X), [T0sat; T1sat; T2sat; 800; 800; 800])
n3_ss=X_initial(6)
L4_ss=((c3.*sqrt(X_initial(6)))./2)
L0_ss=(c1.*sqrt(X_initial(4)))-(alpha.*(X_initial(1)-T0sat))

New values and parameters
K4=0.05;        %1/s
K0=0.5*10^6;    %J/mol
% L4=L4_ss+K4*(Y(6)-n3_ss)  Case A
% Q_new = Q_ss +K0*(((c1.*sqrt(Y1(4))-alpha.*(Y1(1)-T0sat)))-L0_ss)    Case A+B
% L3= ((c3.*sqrt(Y(6)))-(L4_ss+K4*(Y(6)-n3_ss)))
% L4=L4_ss                Case B
External variables
LF_change= @(t) 20 + 2 *(t>=60) -2*(t>=600); % mol/s


Non-linear model
% Y(1) = T0
% Y(2) = T1
% Y(3) = T2
% Y(4) = n1
% Y(5) = n2
% Y(6) = n3
% Case A 
dY1_dt = @(t, Y1) [((((c1.*sqrt(Y1(4))).*Y1(2))-((c1.*sqrt(Y1(4)) - alpha.*(Y1(1)-T0sat)).*Y1(1))-((alpha.*(Y1(1)-T0sat)).*(Y1(1)+(gammav/cp)))+((Q_ss+K0*(((c1.*sqrt(Y1(4)))-(alpha.*(Y1(2)-T1sat)))-L0_ss))/cp))/900); ...
                 (c2.*sqrt(Y1(5)))-(c1.*sqrt(Y1(4)))-(alpha.*(Y1(2)-T1sat)); ...
                 (((c2.*sqrt(Y1(5)).*(Y1(3)-Y1(2)))-((alpha.*(Y1(2)-T1sat)).*(gammav/cp))))./(((c1*sqrt(Y1(4)))/c1)^2); ...
                 LF_change(t)+((c3.*sqrt(Y1(6)))-(L4_ss+K4*(Y1(6)-n3_ss)))-(c2.*sqrt(Y1(5)))+(alpha.*(Y1(1)-T0sat))+(alpha.*(Y1(2)-T1sat))-(alpha.*(Y1(3)-T2sat)); ...
                 ((LF_change(t).*(TF-Y1(3)))+((c3.*sqrt(Y1(6)))-(L4_ss+K4*(Y1(6)-n3_ss))).*(T3-Y1(3))+((alpha.*(Y1(1)-T0sat)).*(Y1(1)-Y1(3)+(gammav/cp)))+((alpha.*(Y1(2)-T1sat)).*(Y1(2)-Y1(3)+(gammav/cp)))-((alpha.*(Y1(3)-T2sat)).*(gammav/cp)))./((((c2*sqrt(Y1(5)))/c2)^2)); ...
                 (alpha.*(Y1(3)-T2sat))-((c3.*sqrt(Y1(6)))-(L4_ss+K4*(Y1(6)-n3_ss)))-(L4_ss+K4*(Y1(6)-n3_ss))];

% Case B 
dY2_dt = @(t, Y2) [((((c1.*sqrt(Y2(4))).*Y2(2))-((c1.*sqrt(Y2(4)) - alpha.*(Y2(1)-T0sat)).*Y2(1))-((alpha.*(Y2(1)-T0sat)).*(Y2(1)+(gammav/cp)))+((Q_ss+K0*(((c1.*sqrt(Y2(4)))-(alpha.*(Y2(2)-T1sat)))-L0_ss))/cp))/900); ...
                 (c2.*sqrt(Y2(5)))-(c1.*sqrt(Y2(4)))-(alpha.*(Y2(2)-T1sat)); ...
                 (((c2.*sqrt(Y2(5)).*(Y2(3)-Y2(2)))-((alpha.*(Y2(2)-T1sat)).*(gammav/cp))))./(((c1*sqrt(Y2(4)))/c1)^2); ...
                 LF_change(t)+((c3.*sqrt(Y2(6)))-L4_ss)-(c2.*sqrt(Y2(5)))+(alpha.*(Y2(1)-T0sat))+(alpha.*(Y2(2)-T1sat))-(alpha.*(Y2(3)-T2sat)); ...
                 ((LF_change(t).*(TF-Y2(3)))+((c3.*sqrt(Y2(6)))-L4_ss).*(T3-Y2(3))+((alpha.*(Y2(1)-T0sat)).*(Y2(1)-Y2(3)+(gammav/cp)))+((alpha.*(Y2(2)-T1sat)).*(Y2(2)-Y2(3)+(gammav/cp)))-((alpha.*(Y2(3)-T2sat)).*(gammav/cp)))./((((c2*sqrt(Y2(5)))/c2)^2)); ...
                 (alpha.*(Y2(3)-T2sat))-((c3.*sqrt(Y2(6)))-L4_ss)-L4_ss];



Solve system dynamics
t=linspace(0,3600,10000);
[t,Y1]=ode23s(dY1_dt,t,X_initial);
[t,Y2]=ode23s(dY2_dt,t,X_initial);
%  Case A 
T0_A1=Y1(:,1); 
T1_A1=Y1(:,2);
T2_A1=Y1(:,3);
n1_A1=Y1(:,4);
n2_A1=Y1(:,5);
n3_A1=Y1(:,6);
fig1 = figure(1);
%plot(t,T0_A1,t,T1_A1,t,T2_A1);
%plot(t,n1_A1,t,n2_A1,t,n3_A1);
%  Case B 
T0_A2=Y2(:,1); 
T1_A2=Y2(:,2);
T2_A2=Y2(:,3);
n1_A2=Y2(:,4);
n2_A2=Y2(:,5);
n3_A2=Y2(:,6);
fig3 = figure(3);
%plot(t,T0_A2,t,T1_A2,t,T2_A2);
%plot(t,n1_A2,t,n2_A2,t,n3_A2);