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.
