How to solve a system of ODEs using ode23s with multiple inputs.

The code works if I use one input (input1) where it calls for the numerical integration ([t,u] = ode23s(@(t,u) gene(t,u,k,input1), tspan, init); and if I define all the inputs in the system as in1(t). This occurs in eqs(6),(16) and (37).

Is it possible though, to have 3 different intput instances, so that the second term in eq(6) will be k26*in2(t), the first term of eq(16) will be k(47)*in3(t)*u(21) and the first term in eq(37) will be k(54)*in4(t)*u(35)?

The code used is given below:

clear
clc
% Set the initial values
u(1) = 5;
u(2) = 1;
u(3) = 1;
u(4) = 1;
u(5) = 1;
u(6) = 50;
u(7) = 1;
u(8) = 100;
u(9) = 5;
u(10) = 1;
u(11) = 1;
u(12) = 1;
u(13) = 1;
u(14) = 1;
u(15) = 1;
u(16) = 1;
u(17) = 1;
u(18) = 1;
u(19) = 1;
u(20) = 1;
u(21) = 1;
u(22) = 1;
u(23) = 1;
u(24) = 1;
u(25) = 1;
u(26) = 1;
u(27) = 1;
u(28) = 1;
u(29) = 1;
u(30) = 1;
u(31) = 1;
u(32) = 1;
u(33) = 1;
u(34) = 1;
u(35) = 1;
u(36) = 1;
u(37) = 1;
u(38) = 1;
u(39) = 1;
u(40) = 1;
u(41) = 1;
%% Define anonymous function input1(t)
%  Time: 0<=t<50 50<=t<100  100<=t<150  150<=t<200  200<=t  
%  Input:  0.5      1          1.5          1         0.5
input1 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
input2 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
input3 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
input4 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
% Set the model parameters
k1 = 3;
k2 = 2;
k3 = 5;
k4 = 4;
k5 = 5;
k6 = 1;
k7 = 4;
k8 = 3;
k9 = 1;
k10 = 1.2;
k11 = 1;
k12 = 9;
k13 = 1;
k14 = 1;
k15 = 1;
k16 = 30;
k17 = 30;
k18 = 1.5;
k19 = 25;
k20 = 1;
k21 = 10;
k22 = 20;
k23 = 1;
k24 = 0.5;
k25 = 1;
k26 = 1;
k27 = 1;
k28 = 1;
k29 = 1;
k30 = 1;
k31 = 1;
k32 = 1;
k33 = 1;
k34 = 1;
k35 = 1;
k36 = 1;
k37 = 1;
k38 = 1;
k39 = 1;
k40 = 1;
k41 = 1;
k42 = 1;
k43 = 1;
k44 = 1;
k45 = 1;
k46 = 1;
k47 = 1;
k48 = 1;
k49 = 1;
k50 = 1;
k51 = 1;
k52 = 1;
k53 = 20;
k54 = 1;
k55 = 1;
k56 = 1;
k57 = 1;
k58 = 1;
k59 = 1;
k60 = 1;
k61 = 1;
k62 = 1;
k63 = 1;
k64 = 1;
k65 = 1;
k66 = 1;
k67 = 1;
k68 = 1;
k69 = 1;
k70 = 1;
%% Perform the numerical integration
k = [k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25,k26,k27,k28,k29,k30,k31,k32,k33,k34,k35,k36,k37,k38,k39,k40,k41,k42,k43,k44,k45,k46,k47,k48,k49,k50,k51,k52,k53,k54,k55,k56,k57,k58,k59,k60,k61,k62,k63,k64,k65,k66,k67,k68,k69,k70];
init = [u(1) u(2) u(3) u(4) u(5) u(6) u(7) u(8) u(9) u(10) u(11) u(12) u(13) u(14) u(15) u(16) u(17) u(18) u(19) u(20) u(21) u(22) u(23) u(24) u(25) u(26) u(27) u(28) u(29) u(30) u(31) u(32) u(33) u(34) u(35) u(36) u(37) u(38) u(39) u(40) u(41)];
tspan = [0 250];
[t,u] = ode23s(@(t,u) gene(t,u,k,input1), tspan, init);
%% Plot results
t1=0:.5:250;
tiledlayout(4,1)
nexttile 
hold on
plot(t1,input1(t1),'-k', 'LineWidth',2.0)
title(''); legend('CL(Input)')
xlabel('Time t'); ylabel('Concentration')
hold off
ylim([0 2])
nexttile 
hold on
plot(t,u(:,6),'-',t,u(:,7),'-',t,u(:,8),'-',t,u(:,9),'-',t,u(:,18),'-',t,u(:,39),'-',t,u(:,29),'-', 'LineWidth',2.0)
title(''); legend('Cf','Cp','Ce','E','Ce1','Ce2','B')
xlabel('Time t'); ylabel('Concentration')
ylim([0 4])
nexttile 
hold on
plot(t,u(:,2),'-',t,u(:,5),'-',t,u(:,11),'-', 'LineWidth',2.0)
title(''); legend('C','R','H')
xlabel('Time t'); ylabel('Concentration')
ylim([0 60])
nexttile 
hold on
plot(t,u(:,1),'-',t,u(:,3),'-',t,u(:,10),'-',t,u(:,12),'-',t,u(:,13),'-',t,u(:,4),'-', 'LineWidth',2.0)
title('','FontSize',16); legend('Sci','Sr','HR','Sp','P','Sh')
xlabel('Time t'); ylabel('Concentration')
ylim([0 1.5])
%% ODE System
function eqns = gene(t,u,k,in1)
% Using u = [Sci C Sr Sh R Cf Cp Ce  E  HR  H  Sp  P  Sb  Rb Cf1  Sci1  C1  Sb1  Sh1  Rb1  HR1 Cp1 Ce1  E1  H1  CO  Cg  B  Sci2  C2  Sr2 Sn2  Sp2  R2  P2  Cf2  Cp2  Ce2  H2  HR2] 
% Equation    1  2  3  4 5  6  7  8  9  10  11 12  13 14  15  16   17  18  
% 
%Using k17=p1, k18=p2,k19=p3, k20=mu, k21=eta, k22=alpha, k23=theta
eqns = zeros(41,1); % To start with we have 41 empty equations
eqns(1) = k(60) - k(61)*u(1)*u(2);
eqns(2) = k(70)*u(8) - k(61)*u(1)*u(2);
eqns(3) = k(1)*u(1) - k(2)*u(3);
eqns(4) = k(1)*u(1) - k(10)*u(4);
eqns(5) = k(17)*u(3) - k(11)*u(5) - k(14)*u(13)*u(5);
eqns(6) = k(3)*in1(t) + k(26)*in1(t) - k(4)*u(6);
eqns(7) = k(4)*u(6) - k(5)*u(7) + k(6)*u(8) - k(12)*u(7)*u(27);
eqns(8) = k(5)*u(7) - k(6)*u(8) + k(9)*u(10)*u(11) - k(7)*u(8) + k(8)*u(9) - k(21)*u(8) - k(22)*u(8)*u(28);
eqns(9) = k(7)*u(8) - k(8)*u(9);
eqns(10) = k(18)*u(4) - k(15)*u(10)*u(8);
eqns(11) = k(62) - k(9)*u(10)*u(11);
eqns(12) = k(1)*u(1) - k(13)*u(12);
eqns(13) = k(19)*u(12) - k(16)*u(13);
eqns(14) = k(1)*u(1) - k(27)*u(14);
eqns(15) = k(28)*u(14) - k(29)*u(15);
eqns(16) = k(47)*in1(t)*u(21) - k(30)*u(16);
eqns(17) = k(63) - k(67)*u(17)*u(18);
eqns(18) = k(64)*u(24) - k(67)*u(17)*u(18);
eqns(19) = k(20)*u(17) - k(23)*u(19);
eqns(20) = k(20)*u(17) - k(25)*u(20);
eqns(21) = k(35)*u(19) - k(36)*u(21);
eqns(22) = k(33)*u(20) - k(38)*u(22)*u(24);
eqns(23) = k(30)*u(16) - k(31)*u(23) + k(32)*u(24) - k(37)*u(23);
eqns(24) = k(31)*u(23) - k(32)*u(24) - k(41)*u(24) + k(42)*u(25) + k(34)*u(22)*u(26) - k(40)*u(24);
eqns(25) = k(41)*u(24) - k(42)*u(25);
eqns(26) = k(65) - k(34)*u(22)*u(26);
eqns(27) = k(37)*u(23) - k(12)*u(7)*u(27);
eqns(28) = k(12)*u(7)*u(27) - k(22)*u(8)*u(28);
eqns(29) = k(22)*u(8)*u(28) - k(24)*u(29);
eqns(30) = k(66) - k(68)*u(30)*u(31);
eqns(31) = k(66)*u(39) - k(68)*u(30)*u(31);
eqns(32) = k(48)*u(30) - k(52)*u(32);
eqns(33) = k(48)*u(30) - k(43)*u(33);
eqns(34) = k(48)*u(30) - k(51)*u(34);
eqns(35) = k(53)*u(32) - k(55)*u(35) - k(59)*u(36)*u(35);
eqns(36) = k(57)*u(34) - k(58)*u(36);
eqns(37) = k(54)*in1(t)*u(35) - k(56)*u(37);
eqns(38) = k(56)*u(37) - k(49)*u(38) + k(50)*u(39);
eqns(39) = k(49)*u(38) - k(50)*u(39) + k(46)*u(41)*u(40) - k(39)*u(39);
eqns(40) = k(69) - k(46)*u(41)*u(40);
eqns(41) = k(44)*u(33) - k(45)*u(41)*u(39);
end

[Walter Roberson]

I recommend that you set up the equations using Symbolic Toolbox, and then follow the flow shown in the first example of odeFunction() to create an anonymous function that is suitable for passing to ode23s()

[Ron_S]

Thank you @Walter Roberson, I'm slowly getting there and I think it is going to work. I did a toy system run first and it worked.

Answer 1

[Sam Chak]

Hi @Ron_S

I simply wanted to ensure that the code was functioning correctly with the given external inputs. Therefore, I didn't make any touch-ups or beautify the display. My main focus was on getting the code to run accurately. Perhaps I'll do some brief explanations afterwards.

% Set the initial values
u(1) = 5;
u(2) = 1;
u(3) = 1;
u(4) = 1;
u(5) = 1;
u(6) = 50;
u(7) = 1;
u(8) = 100;
u(9) = 5;
u(10) = 1;
u(11) = 1;
u(12) = 1;
u(13) = 1;
u(14) = 1;
u(15) = 1;
u(16) = 1;
u(17) = 1;
u(18) = 1;
u(19) = 1;
u(20) = 1;
u(21) = 1;
u(22) = 1;
u(23) = 1;
u(24) = 1;
u(25) = 1;
u(26) = 1;
u(27) = 1;
u(28) = 1;
u(29) = 1;
u(30) = 1;
u(31) = 1;
u(32) = 1;
u(33) = 1;
u(34) = 1;
u(35) = 1;
u(36) = 1;
u(37) = 1;
u(38) = 1;
u(39) = 1;
u(40) = 1;
u(41) = 1;

init = [u(1) u(2) u(3) u(4) u(5) u(6) u(7) u(8) u(9) u(10) u(11) u(12) u(13) u(14) u(15) u(16) u(17) u(18) u(19) u(20) u(21) u(22) u(23) u(24) u(25) u(26) u(27) u(28) u(29) u(30) u(31) u(32) u(33) u(34) u(35) u(36) u(37) u(38) u(39) u(40) u(41)];
tspan = [0 250];
[t,u] = ode23s(@(t,u) gene(t,u), tspan, init);

%% Plot results
t1 = 0:.5:250;
input1 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
tiledlayout(4,1)
nexttile 
hold on
plot(t1,input1(t1),'-k', 'LineWidth',2.0)
title(''); legend('CL(Input)')
xlabel('Time t'); ylabel('Concentration')
hold off
ylim([0 2])
nexttile 
hold on
plot(t,u(:,6),'-',t,u(:,7),'-',t,u(:,8),'-',t,u(:,9),'-',t,u(:,18),'-',t,u(:,39),'-',t,u(:,29),'-', 'LineWidth',2.0)
title(''); legend('Cf','Cp','Ce','E','Ce1','Ce2','B')
xlabel('Time t'); ylabel('Concentration')
ylim([0 4])
nexttile 
hold on
plot(t,u(:,2),'-',t,u(:,5),'-',t,u(:,11),'-', 'LineWidth',2.0)
title(''); legend('C','R','H')
xlabel('Time t'); ylabel('Concentration')
ylim([0 60])
nexttile 
hold on
plot(t,u(:,1),'-',t,u(:,3),'-',t,u(:,10),'-',t,u(:,12),'-',t,u(:,13),'-',t,u(:,4),'-', 'LineWidth',2.0)
title('','FontSize',16); legend('Sci','Sr','HR','Sp','P','Sh')
xlabel('Time t'); ylabel('Concentration')
ylim([0 1.5])

%% ODE System
function eqns = gene(t, u)
    % Using u = [Sci C Sr Sh R Cf Cp Ce  E  HR  H  Sp  P  Sb  Rb Cf1  Sci1  C1  Sb1  Sh1  Rb1  HR1 Cp1 Ce1  E1  H1  CO  Cg  B  Sci2  C2  Sr2 Sn2  Sp2  R2  P2  Cf2  Cp2  Ce2  H2  HR2] 
    % Equation    1  2  3  4 5  6  7  8  9  10  11 12  13 14  15  16   17  18  
    % 
    % Using k17=p1, k18=p2,k19=p3, k20=mu, k21=eta, k22=alpha, k23=theta
    
    % Set the model parameters
    k1 = 3;
    k2 = 2;
    k3 = 5;
    k4 = 4;
    k5 = 5;
    k6 = 1;
    k7 = 4;
    k8 = 3;
    k9 = 1;
    k10 = 1.2;
    k11 = 1;
    k12 = 9;
    k13 = 1;
    k14 = 1;
    k15 = 1;
    k16 = 30;
    k17 = 30;
    k18 = 1.5;
    k19 = 25;
    k20 = 1;
    k21 = 10;
    k22 = 20;
    k23 = 1;
    k24 = 0.5;
    k25 = 1;
    k26 = 1;
    k27 = 1;
    k28 = 1;
    k29 = 1;
    k30 = 1;
    k31 = 1;
    k32 = 1;
    k33 = 1;
    k34 = 1;
    k35 = 1;
    k36 = 1;
    k37 = 1;
    k38 = 1;
    k39 = 1;
    k40 = 1;
    k41 = 1;
    k42 = 1;
    k43 = 1;
    k44 = 1;
    k45 = 1;
    k46 = 1;
    k47 = 1;
    k48 = 1;
    k49 = 1;
    k50 = 1;
    k51 = 1;
    k52 = 1;
    k53 = 20;
    k54 = 1;
    k55 = 1;
    k56 = 1;
    k57 = 1;
    k58 = 1;
    k59 = 1;
    k60 = 1;
    k61 = 1;
    k62 = 1;
    k63 = 1;
    k64 = 1;
    k65 = 1;
    k66 = 1;
    k67 = 1;
    k68 = 1;
    k69 = 1;
    k70 = 1;
    %% Perform the numerical integration
    k = [k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25,k26,k27,k28,k29,k30,k31,k32,k33,k34,k35,k36,k37,k38,k39,k40,k41,k42,k43,k44,k45,k46,k47,k48,k49,k50,k51,k52,k53,k54,k55,k56,k57,k58,k59,k60,k61,k62,k63,k64,k65,k66,k67,k68,k69,k70];
    
    %% Define anonymous function input1(t)
    %  Time: 0<=t<50 50<=t<100  100<=t<150  150<=t<200  200<=t  
    %  Input:  0.5      1          1.5          1         0.5
    in1 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
    in2 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
    in3 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
    in4 = @(t) 0.5 + 0.5*(50<=t & t<100) + 1*(100<=t & t<150) + 0.5*(150<=t & t<200);
    
    eqns     = zeros(41,1); % To start with we have 41 empty equations
    eqns(1)  = k(60) - k(61)*u(1)*u(2);
    eqns(2)  = k(70)*u(8) - k(61)*u(1)*u(2);
    eqns(3)  = k(1)*u(1) - k(2)*u(3);
    eqns(4)  = k(1)*u(1) - k(10)*u(4);
    eqns(5)  = k(17)*u(3) - k(11)*u(5) - k(14)*u(13)*u(5);
    eqns(6)  = k(3)*in1(t) + k(26)*in2(t) - k(4)*u(6);
    eqns(7)  = k(4)*u(6) - k(5)*u(7) + k(6)*u(8) - k(12)*u(7)*u(27);
    eqns(8)  = k(5)*u(7) - k(6)*u(8) + k(9)*u(10)*u(11) - k(7)*u(8) + k(8)*u(9) - k(21)*u(8) - k(22)*u(8)*u(28);
    eqns(9)  = k(7)*u(8) - k(8)*u(9);
    eqns(10) = k(18)*u(4) - k(15)*u(10)*u(8);
    eqns(11) = k(62) - k(9)*u(10)*u(11);
    eqns(12) = k(1)*u(1) - k(13)*u(12);
    eqns(13) = k(19)*u(12) - k(16)*u(13);
    eqns(14) = k(1)*u(1) - k(27)*u(14);
    eqns(15) = k(28)*u(14) - k(29)*u(15);
    eqns(16) = k(47)*in3(t)*u(21) - k(30)*u(16);
    eqns(17) = k(63) - k(67)*u(17)*u(18);
    eqns(18) = k(64)*u(24) - k(67)*u(17)*u(18);
    eqns(19) = k(20)*u(17) - k(23)*u(19);
    eqns(20) = k(20)*u(17) - k(25)*u(20);
    eqns(21) = k(35)*u(19) - k(36)*u(21);
    eqns(22) = k(33)*u(20) - k(38)*u(22)*u(24);
    eqns(23) = k(30)*u(16) - k(31)*u(23) + k(32)*u(24) - k(37)*u(23);
    eqns(24) = k(31)*u(23) - k(32)*u(24) - k(41)*u(24) + k(42)*u(25) + k(34)*u(22)*u(26) - k(40)*u(24);
    eqns(25) = k(41)*u(24) - k(42)*u(25);
    eqns(26) = k(65) - k(34)*u(22)*u(26);
    eqns(27) = k(37)*u(23) - k(12)*u(7)*u(27);
    eqns(28) = k(12)*u(7)*u(27) - k(22)*u(8)*u(28);
    eqns(29) = k(22)*u(8)*u(28) - k(24)*u(29);
    eqns(30) = k(66) - k(68)*u(30)*u(31);
    eqns(31) = k(66)*u(39) - k(68)*u(30)*u(31);
    eqns(32) = k(48)*u(30) - k(52)*u(32);
    eqns(33) = k(48)*u(30) - k(43)*u(33);
    eqns(34) = k(48)*u(30) - k(51)*u(34);
    eqns(35) = k(53)*u(32) - k(55)*u(35) - k(59)*u(36)*u(35);
    eqns(36) = k(57)*u(34) - k(58)*u(36);
    eqns(37) = k(54)*in4(t)*u(35) - k(56)*u(37);
    eqns(38) = k(56)*u(37) - k(49)*u(38) + k(50)*u(39);
    eqns(39) = k(49)*u(38) - k(50)*u(39) + k(46)*u(41)*u(40) - k(39)*u(39);
    eqns(40) = k(69) - k(46)*u(41)*u(40);
    eqns(41) = k(44)*u(33) - k(45)*u(41)*u(39);
end

[Ron_S]

Thank you so much @Sam Chak, this is exactly what I wanted to see. I'm learning how to do this in a more succinct and efficient way as suggested by @Walter Roberson, and I think I'm slowly getting there. Thank you to both contributors.