Runge-Kutta for multiple variables

Hello everyone, I am developing a code for resolving a kinetic model of reactions and calculate the concentrations of 4 compounds. I'm implementing the method 4th order Runge-Kutta for multiple variables to calculate each concentration. For this I am coding 4 separate loops for each concentration. I am having trouble with the response vectors; it seems that the course of the reaction isn't correct. I'll appreciate any help solving this issue also an explanation why this is happening will be really helpful to improve my knowledge in Matlab.

clear all; 
close all; 
clc; 
k = [0.75,0.71,0.68]; 
E0 = 0.5; %Enzima
S0 = 1; %Ssutrato 
C0 = 0;%Complejo 
P0 = 0; %producto 
dt = 10e-4; 
x = 0:dt:3 ; %minutos

E = zeros(1,length(x));
S = zeros(1,length(x));
C = zeros(1,length(x));
P = zeros(1,length(x)); 

fE = @(t,E,S,C)-(k(1).*E.*S)+((k(2)+k(3)).*C);
fS = @(t,E,S,C)-(k(1).*S.*E)+((k(2).*C));
fC = @(t,E,S,C)(k(1).*S.*E)-((k(2)+k(3)).*C);
fP = @(t,C)k(3).*C; 
 
RUNKE-KUTTA 4 
%Enzima

E(1) = E0; 
for i=1:(length(x)-1)                            
    k_1 = fE(x(i),E(i),S(i),C(i));
    k_2 = fE(x(i)+0.5*dt,E(i)+0.5*dt*k_1,S(i),C(i));
    k_3 = fE((x(i)+0.5*dt),(E(i)+0.5*dt*k_2),(S(i)),C(i));
    k_4 = fE((x(i)+dt),(E(i)+k_3*dt),(S(i)),(C(i)));
    E(i+1) = E(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt;
end
% 
% %Sustrato 

S(1) = S0; 
for i=1:(length(x)-1)                            
    k_1 = fS(x(i),E(i),S(i),C(i));
    k_2 = fS(x(i)+0.5*dt,E(i),S(i)+0.5*dt*k_1,C(i));
    k_3 = fS((x(i)+0.5*dt),E(i),(S(i)+0.5*dt*k_2),C(i));
    k_4 = fS((x(i)+dt),E(i),(S(i)+k_3*dt),C(i));
    S(i+1) = S(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt; 
end
% 
% %Complejo 

C(1) = C0;
for i=1:(length(x)-1)                            
    k_1 = fC(x(i),C(i),S(i),E(i));
    k_2 = fC(x(i)+0.5*dt,E(i),S(i),C(i)+0.5*dt*k_1);
    k_3 = fC((x(i)+0.5*dt),E(i),S(i),(C(i)+0.5*dt*k_2));
    k_4 = fC((x(i)+dt),E(i),S(i),(C(i)+k_3*dt));
    C(i+1) = C(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt; 
end
% 
% %Producto 

P(1) = P0; 
for i=1:(length(x)-1)                            
    k_1 = fP(x(i),P(i));
    k_2 = fP(x(i)+0.5*dt,P(i)+0.5*dt*k_1);
    k_3 = fP((x(i)+0.5*dt),(P(i)+0.5*dt*k_2));
    k_4 = fP((x(i)+dt),(P(i)+k_3*dt));
    P(i+1) = P(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt; 
end

figure(1)
plot(x,E,x,S,x,C,x,P)
xlabel('$t (min)$',Interpreter='latex')
ylabel('Concentraciones $(\frac{mmol}{L})$',Interpreter='latex')
legend('Enzima','Sustrato','Compuesto','Producto')
grid on 
grid minor 

Answer 1

[Davide Masiello]

Hi @Juliana Quintana, the problem is that by defining 4 different function handles and 4 different loops, you have decoupled the ode system.

You can fix this as shown below

clear,clc

k   = [0.75,0.71,0.68]; 
E0  = 0.5;                  % Enzima
S0  = 1;                    % Sustrato 
C0  = 0;                    % Complejo 
P0  = 0;                    % Producto 
h   = 1e-4; 

t = 0:h:3 ;                 % Minutos

y       = zeros(4,length(t));
y(:,1)  = [E0;S0;C0;P0];
 
% RUNGE-KUTTA 4
for i = 1:length(t)-1                          
    k1 = odeSystem(t(i),y(:,i),k);
    k2 = odeSystem(t(i)+h/2,y(:,i)+h*k1/2,k);
    k3 = odeSystem(t(i)+h/2,y(:,i)+h*k2/2,k);
    k4 = odeSystem(t(i)+h,y(:,i)+h*k3,k);
    y(:,i+1) = y(:,i) + (1/6)*(k1+2*k2+2*k3+k4)*h;
end

E = y(1,:);
S = y(2,:);
C = y(3,:);
P = y(4,:);

figure(1)
plot(t,E,t,S,t,C,t,P)
xlabel('$t (min)$',Interpreter='latex')
ylabel('Concentraciones $(\frac{mmol}{L})$',Interpreter='latex')
legend('Enzima','Sustrato','Compuesto','Producto')
grid on 
grid minor

function dydt = odeSystem(t,y,k)
E = y(1);
S = y(2);
C = y(3);
P = y(4);

dydt(1,1) = -(k(1).*E.*S)+((k(2)+k(3)).*C);
dydt(2,1) = -(k(1).*S.*E)+((k(2).*C));
dydt(3,1) = (k(1).*S.*E)-((k(2)+k(3)).*C);
dydt(4,1) = k(3).*C; 
end