MATLAB: Step change for ODE

MATLABodestep change

Hi,
I have a function that solves the dynamic behavior of a reactor over time given certain input conditions using the ODE solver. However, I would like to be able to see what would happen if a step change in the input condition or some other variable is applied at a certain t. I believe Simulink can be used to determine the response to a step function, but I'm not sure how to pass my function to Simulink. For example, in my code below, the input condition is Tc. I would like to see the change if the variable q was changed from 150 to 100 at t = 20 if Tc is the same or if Tc was changed at t = 20.
Thanks for any help!
function ivp
% Inputs
tval=0:0.01:20;
Tc = [290; 300; 305; 310];
% Differential Equation
Ca0=0.5; T0=350; Tc0 = 300;
Int = [Ca0; T0];
[t,Y1] = ode15s(@ode,tval,Int,[],Tc(1));
% Plots
figure()
plot(t,Y1(:,2),'r')
xlabel('Time (min)'); ylabel('Reactor temperature (K)');legend('290K');
figure()
plot(t,Y1(:,1),'r')
xlabel('Time (min)');ylabel('Reactant A concentration (mol/L)');legend('290K');
end
function [Dy] = ode(t,y,Tc)
q=150; Cai=1; Ti=350; V=100; p=1000; C=0.239; H=5*10^4;
ER=8750; k0=7.2*10^10; UA=5*10^4; 
% y(1) is Ca   Dy(1) is deriv
% y(2) is T    Dy(2) is deriv
Dy  = zeros(2,1);
Dy(1) = q/V*(Cai-y(1))-k0*exp(-ER/y(2))*y(1);
Dy(2) = q/V*(Ti-y(2))+H/(p*C)*k0*exp(-ER/y(2))*y(1)+UA/(V*p*C)*(Tc-y(2));
end

Best Answer

Try this: 
function ivp
% Inputs

tval=0:0.01:20;
Tc = [290; 300; 305; 310];
% Differential Equation

Ca0=0.5; T0=350; Tc0 = 300;
Int = [Ca0; T0];
[t1,Y11] = ode15s(@ode,tval,Int,[],Tc(1));
tval = tval+20.01;                                              % New ‘t’

Int = Y11(end,:);                                               % Last Values Are New Initial Conditions

[t2,Y12] = ode15s(@ode,tval,Int,[],Tc(2));                      % New Value For ‘Tc’

t = [t1; t2];                                                   % Concatenate Time Vectors
Y1 = [Y11; Y12];                                                % Concatenate Integrated Outputs
% Plots

figure()
plot(t,Y1(:,2),'r')
xlabel('Time (min)'); ylabel('Reactor temperature (K)');legend('290K');
figure()
plot(t,Y1(:,1),'r')
xlabel('Time (min)');ylabel('Reactant A concentration (mol/L)');legend('290K');
end
function [Dy] = ode(t,y,Tc)
% q=150; 
q = 150.*(t<=20) + 100.*(t>20);                                 % Define ‘q’ With Step Change At ‘t=20’
Cai=1; Ti=350; V=100; p=1000; C=0.239; H=5*10^4;
q = 150.*(t<=20) + 100.*(t>20);
ER=8750; k0=7.2*10^10; UA=5*10^4; 
% y(1) is Ca   Dy(1) is deriv
% y(2) is T    Dy(2) is deriv
Dy  = zeros(2,1);
Dy(1) = q/V*(Cai-y(1))-k0*exp(-ER/y(2))*y(1);
Dy(2) = q/V*(Ti-y(2))+H/(p*C)*k0*exp(-ER/y(2))*y(1)+UA/(V*p*C)*(Tc-y(2));
end
This changes ‘Tc’ in the second ode15s call, and changes ‘q’ inside the ‘ode’ function. Experiment with other options.
EDIT — (9 Aug 2019 at 21:59)
To plot the temperatures and their effects on the plot, with corresponding legend descriptions, the code changes: 
function ivp
% Inputs
tval=0:0.01:20;
Tc = [290; 300; 305; 310];
% Differential Equation
Ca0=0.5; T0=350; Tc0 = 300;
Int = [Ca0; T0];
[t1,Y11] = ode15s(@ode,tval,Int,[],Tc(1));
tval = tval+20.01;                                              % New ‘t’
Int = Y11(end,:);                                               % Last Values Are New Initial Conditions
[t2,Y12] = ode15s(@ode,tval,Int,[],Tc(2));                      % New Value For ‘Tc’
% t = [t1; t2];                                                   % Concatenate Time Vectors
% Y1 = [Y11; Y12];                                                % Concatenate Integrated Outputs
% Plots
figure()
plot(t1,Y11(:,2),'r')
hold on
plot(t2,Y12(:,2),'g')
hold off
xlabel('Time (min)')
ylabel('Reactor temperature (K)')
legend('290K','300K');
figure()
plot(t1,Y11(:,1),'r')
hold on
plot(t2,Y12(:,1),'g')
hold off
xlabel('Time (min)')
ylabel('Reactant A concentration (mol/L)')
legend('290K','300K');
end
If you have more temperatures and more values for ‘q’, a loop will be most efficient, and a different definition for ‘q’ will be necessary. This code works for two temperatures and two values for ‘q’.
Except for the plot section changes, my code is unchanged fro m my original Answer.
Related Question