Vectorizing a Large Function with a Command

Hello, everyone. I am trying to run this program, however, I have run into a big problem. As you can see below, the code is used to solve a 2DOF system of differential equations for the dissociation reaction of CO2. However, the differential equations are much too large to tediously sit through and change each operation in order to allow for solving the system as an array. I am wondering, there has to be a way to automate this process. To automatically "Vectorize" a function so it can be solved with matrices. In engineering, there are countless complex equations that would take considerable effort to change each operation to suite the array solver needs. I have also been running into a strange issue where the command window gives a "Too many input arguments" error for the function of dy1, when the variables chosen for the runge kutta integration for loop are all in the function itself.

I am very very close to solving this reaction, and doing so would help me move further along in my internal combustion engine research. I am open to any suggestions, and any help would be greatly appreciated. Thanks.

P

Code:

% Chemical Equilibrium

% Solve the following chemical equation: 
% aCO2 -> n1[CO2]+n2[O2]+n3[CO]+n4[O]

% Goal: We want to be able to solve the chemical reaction above for all
% products assuming chemical equilibrium

% Assume a fraction of CO2 and O2 dissociates. We will call these parameters
% "y1" and "y2" 

% Write out the mass balance: 
% a -> a(1-y1)[CO2]+(1/2)ay1(1-y2)[O2]+ay1[CO]+ay1[CO]+ay1y2[O]

% From this, we can derive the total number of products:

% np = a[1-y1+(1/2)y1(1-y2)+y1+y1y2]
% np = (1/2)a[2+y1(1+y2)]

% From this, we define each mole fraction

% x1 = 2(1-y1)/(2+y1(1+y2))
% x2 = y1(1-y2)/(2+y1(1+y2))
% x3 = 2y1/(2+y1(1+y2))
% x4 = 2y1y2/(2+y1(1+y2))

% From here, we can plug these values into the equilibrium constant
% equation and solve for the time derivate of y. To solve for the
% differential equation, initialization parameters xij are used.

clearvars

syms K1 dK1 K2 dK2 y1 dy1 y2 dy2 P dP

% Initialization Parameters: These mole fraction and dissociation constant
% equations are being used to solve for the time derivatives of y1 and y2

x1j = 2*(1-y1)/(2+y1*(1+y2));

x2j = y1*(1-y2)/(2+y1*(1+y2));

x3j = y1*(1-y2)/(2+y1*(1+y2));

x4j = y1*(1-y2)/(2+y1*(1+y2));

% Differential Initialization Parameters

dx1j = diff(x1j,y1)*dy1+diff(x1j,y2)*dy2;

dx2j = diff(x2j,y1)*dy1+diff(x2j,y2)*dy2;

dx3j = diff(x3j,y1)*dy1+diff(x3j,y2)*dy2;

dx4j = diff(x4j,y1)*dy1+diff(x4j,y2)*dy2;

% Dissociation Constant Initialization Parameters

x2jk = (x4j^2*P)/(K2^2);

x3jk = (x1j^2)/(x2j^2*K1^2*P);

%solving the combined equations for dy1 and dy2

eqn = [diff(x2j,y1)*dy1+diff(x2j,y2)*dy2 == diff(x2jk,y1)*dy1 + diff(x2jk,y2)*dy2 + diff(x2jk,P)*dP + diff(x2jk,K1)*dK1+diff(x2jk,K2)*dK2, diff(x3j,y1)*dy1+diff(x3j,y2)*dy2 == diff(x3jk,y1)*dy1 + diff(x3jk,y2)*dy2 + diff(x3jk,P)*dP + diff(x3jk,K1)*dK1+diff(x3jk,K2)*dK2];

sol = solve(eqn,[dy1 dy2]);

dy1 = simplify(sol.dy1,"Steps",50);

dy2 = simplify(sol.dy2,"Steps",50);



%Now that we have equations for dy1 and dy2, we can proceed with the
%numerical analysis of the system

clearvars

h = 0.01;
tf = 60;
N = tf/h;

t = linspace(0,tf,N+1);

y1 = zeros(1,N+1);
y2 = zeros(1,N+1);

y1(1) = 0;
y2(1) = 0; 

dK1 = 1;
dK2 = 1;

K1 = 1;
K2 = 1; 

dP = 1;

P = 1;

for n=1:N
    k1y1 = h*f(y1(n),K1,K2,dK1,dK2,dP,P,y1,y2);
    k2y1 = h*f(y1(n)+1/2*k1y1,K1,K2,dK1,dK2,dP,P,y1,y2);
    k3y1 = h*f(y1(n)+1/2*k2y1,K1,K2,dK1,dK2,dP,P,y1,y2);
    k4y1 = h*f(y1(n)+k3y1,K1,K2,dK1,dK2,dP,P,y1,y2);
    y1(n+1) = y1(n)+(1/6)*(k1y1+2*(k2y1+ k3y1)+ k4y1);
end

for n=1:N
    k1y2 = h*g(y2(n),K1,K2,dK1,dK2,dP,P,y1,y2);
    k2y2 = h*g(y2(n)+1/2*k1y2,K1,K2,dK1,dK2,dP,P,y1,y2);
    k3y2 = h*g(y2(n)+1/2*k2y2,K1,K2,dK1,dK2,dP,P,y1,y2);
    k4y2 = h*g(y2(n)+k3y2,K1,K2,dK1,dK2,dP,P,y1,y2);
    y2(n+1) = y2(n)+(1/6)*(k1y2+2*(k2y2+k3y2)+k4y2);
end

nargin('dy1')


function dy1 = f(K1,K2,dK1,dK2,dP,P,y1,y2)

dy1 = (dP*K1^3*K2*P^2*y1^5*y2^4 - 4*dP*K1^3*K2*P^2*y1^5*y2^3 + 6*dP*K1^3*K2*P^2*y1^5*y2^2 - 4*dP*K1^3*K2*P^2*y1^5*y2 + dP*K1^3*K2*P^2*y1^5 + dP*K1^3*K2*P^2*y1^4*y2^4 - 4*dP*K1^3*K2*P^2*y1^4*y2^3 + 6*dP*K1^3*K2*P^2*y1^4*y2^2 - 4*dP*K1^3*K2*P^2*y1^4*y2 + dP*K1^3*K2*P^2*y1^4 - 2*dK2*K1^3*P^3*y1^5*y2^4 + 8*dK2*K1^3*P^3*y1^5*y2^3 - 12*dK2*K1^3*P^3*y1^5*y2^2 + 8*dK2*K1^3*P^3*y1^5*y2 - 2*dK2*K1^3*P^3*y1^5 - 2*dK2*K1^3*P^3*y1^4*y2^4 + 8*dK2*K1^3*P^3*y1^4*y2^3 - 12*dK2*K1^3*P^3*y1^4*y2^2 + 8*dK2*K1^3*P^3*y1^4*y2 - 2*dK2*K1^3*P^3*y1^4 + 4*dP*K1*K2^3*y1^5*y2^2 + 8*dP*K1*K2^3*y1^5*y2 + 4*dP*K1*K2^3*y1^5 - 4*dP*K1*K2^3*y1^4*y2^2 + 8*dP*K1*K2^3*y1^4*y2 + 12*dP*K1*K2^3*y1^4 - 4*dP*K1*K2^3*y1^3*y2^2 - 24*dP*K1*K2^3*y1^3*y2 - 4*dP*K1*K2^3*y1^3 + 4*dP*K1*K2^3*y1^2*y2^2 - 8*dP*K1*K2^3*y1^2*y2 - 28*dP*K1*K2^3*y1^2 + 16*dP*K1*K2^3*y1*y2 + 16*dP*K1*K2^3 - 4*dP*K1*K2*P*y1^5*y2^3 + 4*dP*K1*K2*P*y1^5*y2^2 + 4*dP*K1*K2*P*y1^5*y2 - 4*dP*K1*K2*P*y1^5 + 8*dP*K1*K2*P*y1^4*y2^3 - 16*dP*K1*K2*P*y1^4*y2^2 + 8*dP*K1*K2*P*y1^4*y2 - 4*dP*K1*K2*P*y1^3*y2^3 + 20*dP*K1*K2*P*y1^3*y2^2 - 28*dP*K1*K2*P*y1^3*y2 + 12*dP*K1*K2*P*y1^3 - 8*dP*K1*K2*P*y1^2*y2^2 + 16*dP*K1*K2*P*y1^2*y2 - 8*dP*K1*K2*P*y1^2 + 8*dK2*K1*P^2*y1^5*y2^3 + 8*dK2*K1*P^2*y1^5*y2^2 - 8*dK2*K1*P^2*y1^5*y2 - 8*dK2*K1*P^2*y1^5 - 16*dK2*K1*P^2*y1^4*y2^3 + 16*dK2*K1*P^2*y1^4*y2^2 + 16*dK2*K1*P^2*y1^4*y2 - 16*dK2*K1*P^2*y1^4 + 8*dK2*K1*P^2*y1^3*y2^3 - 56*dK2*K1*P^2*y1^3*y2^2 + 24*dK2*K1*P^2*y1^3*y2 + 24*dK2*K1*P^2*y1^3 + 32*dK2*K1*P^2*y1^2*y2^2 - 64*dK2*K1*P^2*y1^2*y2 + 32*dK2*K1*P^2*y1^2 + 32*dK2*K1*P^2*y1*y2 - 32*dK2*K1*P^2*y1 + 8*dK1*K2^3*P*y1^5*y2^2 + 16*dK1*K2^3*P*y1^5*y2 + 8*dK1*K2^3*P*y1^5 - 8*dK1*K2^3*P*y1^4*y2^2 + 16*dK1*K2^3*P*y1^4*y2 + 24*dK1*K2^3*P*y1^4 - 8*dK1*K2^3*P*y1^3*y2^2 - 48*dK1*K2^3*P*y1^3*y2 - 8*dK1*K2^3*P*y1^3 + 8*dK1*K2^3*P*y1^2*y2^2 - 16*dK1*K2^3*P*y1^2*y2 - 56*dK1*K2^3*P*y1^2 + 32*dK1*K2^3*P*y1*y2 + 32*dK1*K2^3*P + 16*dK1*K2*P^2*y1^5*y2^2 - 16*dK1*K2*P^2*y1^5 - 16*dK1*K2*P^2*y1^4*y2^2 + 32*dK1*K2*P^2*y1^4*y2 - 16*dK1*K2*P^2*y1^4 - 16*dK1*K2*P^2*y1^3*y2^2 - 32*dK1*K2*P^2*y1^3*y2 + 48*dK1*K2*P^2*y1^3 + 16*dK1*K2*P^2*y1^2*y2^2 - 32*dK1*K2*P^2*y1^2*y2 + 16*dK1*K2*P^2*y1^2 + 32*dK1*K2*P^2*y1*y2 - 32*dK1*K2*P^2*y1)/(16*K1*K2*P*(y1 - 1)*(y1 + y1*y2 + 2)*(K2^2*y1 - 2*P*y1 + 2*K2^2 + 2*P*y1*y2 + K2^2*y1*y2));

end 

function dy2 = g(K1,K2,dK1,dK2,dP,P,y1,y2)

dy2 = -((y2 - 1)*(dP*K1^3*K2*P^2*y1^4*y2^4 - 4*dP*K1^3*K2*P^2*y1^4*y2^3 + 6*dP*K1^3*K2*P^2*y1^4*y2^2 - 4*dP*K1^3*K2*P^2*y1^4*y2 + dP*K1^3*K2*P^2*y1^4 - 2*dK2*K1^3*P^3*y1^4*y2^4 + 8*dK2*K1^3*P^3*y1^4*y2^3 - 12*dK2*K1^3*P^3*y1^4*y2^2 + 8*dK2*K1^3*P^3*y1^4*y2 - 2*dK2*K1^3*P^3*y1^4 + 4*dP*K1*K2^3*y1^4*y2^2 + 8*dP*K1*K2^3*y1^4*y2 + 4*dP*K1*K2^3*y1^4 - 8*dP*K1*K2^3*y1^3*y2^2 + 8*dP*K1*K2^3*y1^3 + 4*dP*K1*K2^3*y1^2*y2^2 - 24*dP*K1*K2^3*y1^2*y2 - 12*dP*K1*K2^3*y1^2 + 16*dP*K1*K2^3*y1*y2 - 16*dP*K1*K2^3*y1 + 16*dP*K1*K2^3 + 4*dP*K1*K2*P*y1^4*y2^3 + 12*dP*K1*K2*P*y1^4*y2^2 - 4*dP*K1*K2*P*y1^4*y2 - 12*dP*K1*K2*P*y1^4 - 4*dP*K1*K2*P*y1^3*y2^3 - 4*dP*K1*K2*P*y1^3*y2^2 + 20*dP*K1*K2*P*y1^3*y2 - 12*dP*K1*K2*P*y1^3 - 8*dP*K1*K2*P*y1^2*y2^2 - 16*dP*K1*K2*P*y1^2*y2 + 24*dP*K1*K2*P*y1^2 - 8*dK2*K1*P^2*y1^4*y2^3 - 8*dK2*K1*P^2*y1^4*y2^2 + 8*dK2*K1*P^2*y1^4*y2 + 8*dK2*K1*P^2*y1^4 + 8*dK2*K1*P^2*y1^3*y2^3 - 24*dK2*K1*P^2*y1^3*y2^2 - 8*dK2*K1*P^2*y1^3*y2 + 24*dK2*K1*P^2*y1^3 + 32*dK2*K1*P^2*y1^2*y2^2 - 32*dK2*K1*P^2*y1^2*y2 + 32*dK2*K1*P^2*y1*y2 - 32*dK2*K1*P^2*y1 + 8*dK1*K2^3*P*y1^4*y2^2 + 16*dK1*K2^3*P*y1^4*y2 + 8*dK1*K2^3*P*y1^4 - 16*dK1*K2^3*P*y1^3*y2^2 + 16*dK1*K2^3*P*y1^3 + 8*dK1*K2^3*P*y1^2*y2^2 - 48*dK1*K2^3*P*y1^2*y2 - 24*dK1*K2^3*P*y1^2 + 32*dK1*K2^3*P*y1*y2 - 32*dK1*K2^3*P*y1 + 32*dK1*K2^3*P + 16*dK1*K2*P^2*y1^4*y2^2 - 16*dK1*K2*P^2*y1^4 - 32*dK1*K2*P^2*y1^3*y2^2 + 32*dK1*K2*P^2*y1^3*y2 + 16*dK1*K2*P^2*y1^2*y2^2 - 64*dK1*K2*P^2*y1^2*y2 + 48*dK1*K2*P^2*y1^2 + 32*dK1*K2*P^2*y1*y2 - 32*dK1*K2*P^2*y1))/(16*K1*K2*P*y1*(y1 - 1)*(y1 + y1*y2 + 2)*(K2^2*y1 - 2*P*y1 + 2*K2^2 + 2*P*y1*y2 + K2^2*y1*y2));

end 

[Jan]

@Jason Louison: "where the command window gives a "Too many input arguments" error for the function of dy1" - Do not improve a failing function, but fix the error at first. Whenever you mention an error in the forum, attach a copy of the complete message. Then the readers do not have to guess the details.

For the vectorization, replacing all * by .* and ^ by .^ should be sufficient already.

[Jason Louison]

Is there a way to do this in one step or command? Again, I’m trying so this in the most efficient way, and the differential equations are very long. Something similar to a “substitute” function, but doing it with the operators instead of the variables.

[Torsten]

You can't solve for y1 and y2 separately - the differential equations are not independent from each other.

[Jason Louison]

I am not solving them explicitly. I am solving the differential equations of their time derivates. By numerically integrating with the 4th Order Runge-Kutta method, I can use the rate of change of y1 and y2 to estimate the actual value of y1 and/or y2 at the next time step. I have done this with matlab before, and will attach a working code to demostrate:

% Chemical Equilibrium Reaction Solving: By Jason Louison

% Solve the following chemical equation: aO2 -> n1[O2]+n2[O]

% Goal: We want to be able to solve the chemical reaction above for all
% products assuming chemical equilibrium

% Assume a fraction of O2 dissociates. We will call this parameter "y" 

%Write out the mass balance: a -> a(1-y)[O2]+2ay*phi[O]

%from this, we can derive the total number of products:

%np a(1-y)+2ay = a(1+y) 

%From this, we define each mole fraction

% x1 = (1-y)/(1+y)

% x2 = (2y)/(1+y)

%From here, we can plug these values into the equilibrium constant
%equation and solve for the time derivate of y. To solve for the
%differential equation, initialization parameters xij are used.

%clearvars

syms K dK y dy P dP

x1i = (1-y)/(1+y);

x2i = (2*y)/(1+y);

dx1i = diff(x1i,y)*dy;

dx2i = diff(x2i,y)*dy;

x2iK = (x1i^2*P)/(K^2);

eqn = diff(x2i,y)*dy == diff(x2iK,y)*dy + diff(x2iK,P)*dP + diff(x2iK,K)*dK;

sol = solve(eqn,dy);

dy = sol; 

%Simulation Parameters: Time step, final time, and Number of Cells

h = 0.01;
tf = 60;
N = tf/h;

%Designation of vector space for time and y

t = linspace(0,tf,N+1);

y = zeros(1,N+1);

%Initial Condition for y

y(1) = 0;

%Designation of vector space for dK, dP, P, and K:

dK = zeros(1,N+1); 

dP = zeros(1,N+1); 

P = zeros(1,N+1);

K = zeros(1,N+1);

%This allows us to vary the above parameters as functions of time. The
%parameters are assigned constant values below

dK = 0.34;

dP = 1;

P = 1;

K = 1;

for n=1:N %fourth-order runge kutta method algorithm
    k1 = h*f(y(n),dK,dP,P,K);
    k2 = h*f(y(n)+1/2*k1,dK,dP,P,K);
    k3 = h*f(y(n)+1/2*k2,dK,dP,P,K);
    k4 = h*f(y(n)+k3,dK,dP,P,K);
    y(n+1) = y(n)+(1/6)*(k1+2*(k2+k3)+k4);
end

%Mole Fractions x1 and x2 as functions of dissociation fraction

x1 = (1-y)./(1+y);

x2 = (2.*y)./(1+y);

%Plot x1, x2, and y vs time

figure(1)
plot(t,x1,"r",t,x2,"b")
title('Dissociation of O_2 vs Time')
xlabel('Time (s)')
ylabel('mole fraction (ni/nt)')
legend({'O_2','O'})

%figure(2)
%plot(t,y)
%title('Dissociation Fraction of O_2 vs Time')
%xlabel('Time (s)')
%ylabel('mole fraction (ni/nt)')

%Main Differential Equation of y: This describes the rate of change of y
%with respect to itself, the time derivatives of Dissociation constant K
%and pressure P (dK and dP), as well as the dissociation constant K and
%pressure P:

function dy = f(y,dK,P,dP,K)

dy = (K.*dP - 2.*P.*dK - K.*dP*y + 2.*P.*dK.*y - K.*dP.*y.^2 + K.*dP.*y.^3 + 2.*P.*dK.*y.^2 - 2.*P.*dK.*y.^3)/(2.*K.^3.*y + 2.*K.^3 + 4.*K.*P - 4.*K.*P.*y);

end 

[Torsten]

Yes, but for the two-equation case, you might have noticed that dy1 depends on y2 and dy2 depends on y1. Thus the two equations cannot be solved independently from each other.

[Jason Louison]

Yes, I understand. Does this mean I need to include both differential equations into one function?

[Torsten]

Yes. Imitate the user interface of the MATLAB ode solvers.

function dy = fun(t,y,pars)

Here, y = [y(1);y(2)] , time t and aybe additional parameters pars are inputs and [dy(1);dy(2)] are outputs.

Your Runge-Kutta code should address this functional interface. Thus in your updates

k1 = h*f(y(n),dK,dP,P,K);
k2 = h*f(y(n)+1/2*k1,dK,dP,P,K);
k3 = h*f(y(n)+1/2*k2,dK,dP,P,K);
k4 = h*f(y(n)+k3,dK,dP,P,K);

f,y(:,n) and ki will all be vectors.

[Jason Louison]

I keep getting a "unrecognized variable 'y' error for the updated code below. I understand how this should work, but for some reason the runge kutta algorithm isn't recognizing 'y' as a parameter.

% Chemical Equilibrium

% Solve the following chemical equation: 
% aCO2 -> n1[CO2]+n2[O2]+n3[CO]+n4[O]

% Goal: We want to be able to solve the chemical reaction above for all
% products assuming chemical equilibrium

% Assume a fraction of CO2 and O2 dissociates. We will call these parameters
% "y1" and "y2" 

% Write out the mass balance: 
% a -> a(1-y1)[CO2]+(1/2)ay1(1-y2)[O2]+ay1[CO]+ay1[CO]+ay1y2[O]

% From this, we can derive the total number of products:

% np = a[1-y1+(1/2)y1(1-y2)+y1+y1y2]
% np = (1/2)a[2+y1(1+y2)]

% From this, we define each mole fraction

% x1 = 2(1-y1)/(2+y1(1+y2))
% x2 = y1(1-y2)/(2+y1(1+y2))
% x3 = 2y1/(2+y1(1+y2))
% x4 = 2y1y2/(2+y1(1+y2))

% From here, we can plug these values into the equilibrium constant
% equation and solve for the time derivate of y. To solve for the
% differential equation, initialization parameters xij are used.

clearvars

%syms K1 dK1 K2 dK2 y1 dy1 y2 dy2 P dP

% Initialization Parameters: These mole fraction and dissociation constant
% equations are being used to solve for the time derivatives of y1 and y2

%x1j = 2*(1-y1)/(2+y1*(1+y2));

%x2j = y1*(1-y2)/(2+y1*(1+y2));

%x3j = y1*(1-y2)/(2+y1*(1+y2));

%x4j = y1*(1-y2)/(2+y1*(1+y2));

% Differential Initialization Parameters

%dx1j = diff(x1j,y1)*dy1+diff(x1j,y2)*dy2;

%dx2j = diff(x2j,y1)*dy1+diff(x2j,y2)*dy2;

%dx3j = diff(x3j,y1)*dy1+diff(x3j,y2)*dy2;

%dx4j = diff(x4j,y1)*dy1+diff(x4j,y2)*dy2;

% Dissociation Constant Initialization Parameters

%x2jk = (x4j^2*P)/(K2^2);

%x3jk = (x1j^2)/(x2j^2*K1^2*P);

%solving the combined equations for dy1 and dy2

%eqn = [diff(x2j,y1)*dy1i+diff(x2j,y2)*dy2i == diff(x2jk,y1)*dy1i + diff(x2jk,y2)*dy2i + diff(x2jk,P)*dPi + diff(x2jk,K1)*dK1i+diff(x2jk,K2)*dK2i, diff(x3j,y1)*dy1+diff(x3j,y2)*dy2 == diff(x3jk,y1)*dy1 + diff(x3jk,y2)*dy2 + diff(x3jk,P)*dP + diff(x3jk,K1)*dK1+diff(x3jk,K2)*dK2];

%sol = solve(eqn,[dy1i dy2i]);

%dy1i = simplify(sol.dy1,"Steps",50);

%dy2i = simplify(sol.dy2,"Steps",50);



%Now that we have equations for dy1 and dy2, we can proceed with the
%numerical analysis of the system

clearvars

h = 0.01;
tf = 60;
N = tf/h;

t = linspace(0,tf,N+1);

y1 = zeros(1,N+1);
y2 = zeros(1,N+1);
dK1 = zeros(1,N+1);
dK2 = zeros(1,N+1);
K1 = zeros(1,N+1);
K2 = zeros(1,N+1);
dP = zeros(1,N+1);
P = zeros(1,N+1);

dK1 = 0;

dK2 = 0;

K1 = 0.5;

K2 = 0.7; 

dP = 0;

P = 1;

for n=1:N
    k1 = h*f(y(n),K1,K2,dK1,dK2,P,dP);
    k2 = h*f(y(n)+1/2*k1,K1,K2,dK1,dK2,P,dP);
    k3 = h*f(y(n)+1/2*k2,K1,K2,dK1,dK2,P,dP);
    k4 = h*f(y(n)+k3,K1,K2,dK1,dK2,P,dP);
    y(n+1) = y(n)+(1/6)*(k1+2*(k2+k3)+k4);
end

%for n=1:N
    %k1y2 = h*g(y2(n),K1,K2,dK1,dK2,dP,P,y1,y2);
    %k2y2 = h*g(y2(n)+1/2*k1y2,K1,K2,dK1,dK2,dP,P,y1,y2);
    %k3y2 = h*g(y2(n)+1/2*k2y2,K1,K2,dK1,dK2,dP,P,y1,y2);
    %k4y2 = h*g(y2(n)+k3y2,K1,K2,dK1,dK2,dP,P,y1,y2);
    %y2(n+1) = y2(n)+(1/6)*(k1y2+2*(k2y2+k3y2)+k4y2);
%end

plot(t,y1,t,y2)

function dy = f(K1,K2,dK1,dK2,dP,P,y)

%Initialization (Initial Values of y1 and y2?)

dy = [0; 0];

dy(1) = (dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^4 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^3 + 6.*dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^2 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2) + dP.*K1.^3.*K2.*P.^2.*y(1).^5 + dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^4 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^3 + 6.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^2 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2) + dP.*K1.^3.*K2.*P.^2.*y(1).^4 - 2.*dK2.*K1.^3.*P.^3.*y(1).^5.*y(2).^4 + 8.*dK2.*K1.^3.*P.^3.*y(1).^5.*y(2).^3 - 12.*dK2.*K1.^3.*P.^3.*y(1).^5.*y(2).^2 + 8.*dK2.*K1.^3.*P.^3.*y(1).^5.*y(2) - 2.*dK2.*K1.^3.*P.^3.*y(1).^5 - 2.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^4 + 8.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^3 - 12.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^2 + 8.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2) - 2.*dK2.*K1.^3.*P.^3.*y(1).^4 + 4.*dP.*K1.*K2.^3.*y(1).^5.*y(2).^2 + 8.*dP.*K1.*K2.^3.*y(1).^5.*y(2) + 4.*dP.*K1.*K2.^3.*y(1).^5 - 4.*dP.*K1.*K2.^3.*y(1).^4.*y(2).^2 + 8.*dP.*K1.*K2.^3.*y(1).^4.*y(2) + 12.*dP.*K1.*K2.^3.*y(1).^4 - 4.*dP.*K1.*K2.^3.*y(1).^3.*y(2).^2 - 24.*dP.*K1.*K2.^3.*y(1).^3.*y(2) - 4.*dP.*K1.*K2.^3.*y(1).^3 + 4.*dP.*K1.*K2.^3.*y(1).^2.*y(2).^2 - 8.*dP.*K1.*K2.^3.*y(1).^2.*y(2) - 28.*dP.*K1.*K2.^3.*y(1).^2 + 16.*dP.*K1.*K2.^3.*y(1).*y(2) + 16.*dP.*K1.*K2.^3 - 4.*dP.*K1.*K2.*P.*y(1).^5.*y(2).^3 + 4.*dP.*K1.*K2.*P.*y(1).^5.*y(2).^2 + 4.*dP.*K1.*K2.*P.*y(1).^5.*y(2) - 4.*dP.*K1.*K2.*P.*y(1).^5 + 8.*dP.*K1.*K2.*P.*y(1).^4.*y(2).^3 - 16.*dP.*K1.*K2.*P.*y(1).^4.*y(2).^2 + 8.*dP.*K1.*K2.*P.*y(1).^4.*y(2) - 4.*dP.*K1.*K2.*P.*y(1).^3.*y(2).^3 + 20.*dP.*K1.*K2.*P.*y(1).^3.*y(2).^2 - 28.*dP.*K1.*K2.*P.*y(1).^3.*y(2) + 12.*dP.*K1.*K2.*P.*y(1).^3 - 8.*dP.*K1.*K2.*P.*y(1).^2.*y(2).^2 + 16.*dP.*K1.*K2.*P.*y(1).^2.*y(2) - 8.*dP.*K1.*K2.*P.*y(1).^2 + 8.*dK2.*K1.*P.^2.*y(1).^5.*y(2).^3 + 8.*dK2.*K1.*P.^2.*y(1).^5.*y(2).^2 - 8.*dK2.*K1.*P.^2.*y(1).^5.*y(2) - 8.*dK2.*K1.*P.^2.*y(1).^5 - 16.*dK2.*K1.*P.^2.*y(1).^4.*y(2).^3 + 16.*dK2.*K1.*P.^2.*y(1).^4.*y(2).^2 + 16.*dK2.*K1.*P.^2.*y(1).^4.*y(2) - 16.*dK2.*K1.*P.^2.*y(1).^4 + 8.*dK2.*K1.*P.^2.*y(1).^3.*y(2).^3 - 56.*dK2.*K1.*P.^2.*y(1).^3.*y(2).^2 + 24.*dK2.*K1.*P.^2.*y(1).^3.*y(2) + 24.*dK2.*K1.*P.^2.*y(1).^3 + 32.*dK2.*K1.*P.^2.*y(1).^2.*y(2).^2 - 64.*dK2.*K1.*P.^2.*y(1).^2.*y(2) + 32.*dK2.*K1.*P.^2.*y(1).^2 + 32.*dK2.*K1.*P.^2.*y(1).*y(2) - 32.*dK2.*K1.*P.^2.*y(1) + 8.*dK1.*K2.^3.*P.*y(1).^5.*y(2).^2 + 16.*dK1.*K2.^3.*P.*y(1).^5.*y(2) + 8.*dK1.*K2.^3.*P.*y(1).^5 - 8.*dK1.*K2.^3.*P.*y(1).^4.*y(2).^2 + 16.*dK1.*K2.^3.*P.*y(1).^4.*y(2) + 24.*dK1.*K2.^3.*P.*y(1).^4 - 8.*dK1.*K2.^3.*P.*y(1).^3.*y(2).^2 - 48.*dK1.*K2.^3.*P.*y(1).^3.*y(2) - 8.*dK1.*K2.^3.*P.*y(1).^3 + 8.*dK1.*K2.^3.*P.*y(1).^2.*y(2).^2 - 16.*dK1.*K2.^3.*P.*y(1).^2.*y(2) - 56.*dK1.*K2.^3.*P.*y(1).^2 + 32.*dK1.*K2.^3.*P.*y(1).*y(2) + 32.*dK1.*K2.^3.*P + 16.*dK1.*K2.*P.^2.*y(1).^5.*y(2).^2 - 16.*dK1.*K2.*P.^2.*y(1).^5 - 16.*dK1.*K2.*P.^2.*y(1).^4.*y(2).^2 + 32.*dK1.*K2.*P.^2.*y(1).^4.*y(2) - 16.*dK1.*K2.*P.^2.*y(1).^4 - 16.*dK1.*K2.*P.^2.*y(1).^3.*y(2).^2 - 32.*dK1.*K2.*P.^2.*y(1).^3.*y(2) + 48.*dK1.*K2.*P.^2.*y(1).^3 + 16.*dK1.*K2.*P.^2.*y(1).^2.*y(2).^2 - 32.*dK1.*K2.*P.^2.*y(1).^2.*y(2) + 16.*dK1.*K2.*P.^2.*y(1).^2 + 32.*dK1.*K2.*P.^2.*y(1).*y(2) - 32.*dK1.*K2.*P.^2.*y(1))/(16.*K1.*K2.*P.*(y(1) - 1).*(y(1) + y(1).*y(2) + 2).*(K2.^2.*y(1) - 2.*P.*y(1) + 2.*K2.^2 + 2.*P.*y(1).*y(2) + K2.^2.*y(1).*y(2)));

dy(2) = -((y(2) - 1).*(dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^4 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^3 + 6.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^2 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2) + dP.*K1.^3.*K2.*P.^2.*y(1).^4 - 2.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^4 + 8.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^3 - 12.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^2 + 8.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2) - 2.*dK2.*K1.^3.*P.^3.*y(1).^4 + 4.*dP.*K1.*K2.^3.*y(1).^4.*y(2).^2 + 8.*dP.*K1.*K2.^3.*y(1).^4.*y(2) + 4.*dP.*K1.*K2.^3.*y(1).^4 - 8.*dP.*K1.*K2.^3.*y(1).^3.*y(2).^2 + 8.*dP.*K1.*K2.^3.*y(1).^3 + 4.*dP.*K1.*K2.^3.*y(1).^2.*y(2).^2 - 24.*dP.*K1.*K2.^3.*y(1).^2.*y(2) - 12.*dP.*K1.*K2.^3.*y(1).^2 + 16.*dP.*K1.*K2.^3.*y(1).*y(2) - 16.*dP.*K1.*K2.^3.*y(1) + 16.*dP.*K1.*K2.^3 + 4.*dP.*K1.*K2.*P.*y(1).^4.*y(2).^3 + 12.*dP.*K1.*K2.*P.*y(1).^4.*y(2).^2 - 4.*dP.*K1.*K2.*P.*y(1).^4.*y(2) - 12.*dP.*K1.*K2.*P.*y(1).^4 - 4.*dP.*K1.*K2.*P.*y(1).^3.*y(2).^3 - 4.*dP.*K1.*K2.*P.*y(1).^3.*y(2).^2 + 20.*dP.*K1.*K2.*P.*y(1).^3.*y(2) - 12.*dP.*K1.*K2.*P.*y(1).^3 - 8.*dP.*K1.*K2.*P.*y(1).^2.*y(2).^2 - 16.*dP.*K1.*K2.*P.*y(1).^2.*y(2) + 24.*dP.*K1.*K2.*P.*y(1).^2 - 8.*dK2.*K1.*P.^2.*y(1).^4.*y(2).^3 - 8.*dK2.*K1.*P.^2.*y(1).^4.*y(2).^2 + 8.*dK2.*K1.*P.^2.*y(1).^4.*y(2) + 8.*dK2.*K1.*P.^2.*y(1).^4 + 8.*dK2.*K1.*P.^2.*y(1).^3.*y(2).^3 - 24.*dK2.*K1.*P.^2.*y(1).^3.*y(2).^2 - 8.*dK2.*K1.*P.^2.*y(1).^3.*y(2) + 24.*dK2.*K1.*P.^2.*y(1).^3 + 32.*dK2.*K1.*P.^2.*y(1).^2.*y(2).^2 - 32.*dK2.*K1.*P.^2.*y(1).^2.*y(2) + 32.*dK2.*K1.*P.^2.*y(1).*y(2) - 32.*dK2.*K1.*P.^2.*y(1) + 8.*dK1.*K2.^3.*P.*y(1).^4.*y(2).^2 + 16.*dK1.*K2.^3.*P.*y(1).^4.*y(2) + 8.*dK1.*K2.^3.*P.*y(1).^4 - 16.*dK1.*K2.^3.*P.*y(1).^3.*y(2).^2 + 16.*dK1.*K2.^3.*P.*y(1).^3 + 8.*dK1.*K2.^3.*P.*y(1).^2.*y(2).^2 - 48.*dK1.*K2.^3.*P.*y(1).^2.*y(2) - 24.*dK1.*K2.^3.*P.*y(1).^2 + 32.*dK1.*K2.^3.*P.*y(1).*y(2) - 32.*dK1.*K2.^3.*P.*y(1) + 32.*dK1.*K2.^3.*P + 16.*dK1.*K2.*P.^2.*y(1).^4.*y(2).^2 - 16.*dK1.*K2.*P.^2.*y(1).^4 - 32.*dK1.*K2.*P.^2.*y(1).^3.*y(2).^2 + 32.*dK1.*K2.*P.^2.*y(1).^3.*y(2) + 16.*dK1.*K2.*P.^2.*y(1).^2.*y(2).^2 - 64.*dK1.*K2.*P.^2.*y(1).^2.*y(2) + 48.*dK1.*K2.*P.^2.*y(1).^2 + 32.*dK1.*K2.*P.^2.*y(1).*y(2) - 32.*dK1.*K2.*P.^2.*y(1)))/(16.*K1.*K2.*P.*y(1).*(y(1) - 1).*(y(1) + y(1).*y(2) + 2).*(K2.^2.*y(1) - 2.*P.*y(1) + 2.*K2.^2 + 2.*P.*y(1).*y(2) + K2.^2.*y(1).*y(2)));

end 

[Jason Louison]

After hours of trying, I have failed to get the code working as you said. My best attempt is shown below, but it gives me a " Unable to perform assignment because the left and right sides have a different number of elements." error. Does anyone know what might be causing this? Thanks in advance.

% Chemical Equilibrium

% Solve the following chemical equation: 
% aCO2 -> n1[CO2]+n2[O2]+n3[CO]+n4[O]

% Goal: We want to be able to solve the chemical reaction above for all
% products assuming chemical equilibrium

% Assume a fraction of CO2 and O2 dissociates. We will call these parameters
% "y1" and "y2" 

% Write out the mass balance: 
% a -> a(1-y1)[CO2]+(1/2)ay1(1-y2)[O2]+ay1[CO]+ay1[CO]+ay1y2[O]

% From this, we can derive the total number of products:

% np = a[1-y1+(1/2)y1(1-y2)+y1+y1y2]
% np = (1/2)a[2+y1(1+y2)]

% From this, we define each mole fraction

% x1 = 2(1-y1)/(2+y1(1+y2))
% x2 = y1(1-y2)/(2+y1(1+y2))
% x3 = 2y1/(2+y1(1+y2))
% x4 = 2y1y2/(2+y1(1+y2))

% From here, we can plug these values into the equilibrium constant
% equation and solve for the time derivate of y. To solve for the
% differential equation, initialization parameters xij are used.

clearvars

%syms K1 dK1 K2 dK2 y1 dy1 y2 dy2 P dP

% Initialization Parameters: These mole fraction and dissociation constant
% equations are being used to solve for the time derivatives of y1 and y2

%x1j = 2*(1-y1)/(2+y1*(1+y2));

%x2j = y1*(1-y2)/(2+y1*(1+y2));

%x3j = y1*(1-y2)/(2+y1*(1+y2));

%x4j = y1*(1-y2)/(2+y1*(1+y2));

% Differential Initialization Parameters

%dx1j = diff(x1j,y1)*dy1+diff(x1j,y2)*dy2;

%dx2j = diff(x2j,y1)*dy1+diff(x2j,y2)*dy2;

%dx3j = diff(x3j,y1)*dy1+diff(x3j,y2)*dy2;

%dx4j = diff(x4j,y1)*dy1+diff(x4j,y2)*dy2;

% Dissociation Constant Initialization Parameters

%x2jk = (x4j^2*P)/(K2^2);

%x3jk = (x1j^2)/(x2j^2*K1^2*P);

%solving the combined equations for dy1 and dy2

%eqn = [diff(x2j,y1)*dy1i+diff(x2j,y2)*dy2i == diff(x2jk,y1)*dy1i + diff(x2jk,y2)*dy2i + diff(x2jk,P)*dPi + diff(x2jk,K1)*dK1i+diff(x2jk,K2)*dK2i, diff(x3j,y1)*dy1+diff(x3j,y2)*dy2 == diff(x3jk,y1)*dy1 + diff(x3jk,y2)*dy2 + diff(x3jk,P)*dP + diff(x3jk,K1)*dK1+diff(x3jk,K2)*dK2];

%sol = solve(eqn,[dy1i dy2i]);

%dy1i = simplify(sol.dy1,"Steps",50);

%dy2i = simplify(sol.dy2,"Steps",50);



%Now that we have equations for dy1 and dy2, we can proceed with the
%numerical analysis of the system

clearvars

h = 0.01;
tf = 60;
N = tf/h;

t = [linspace(0,tf,N+1)]';

dK1 = zeros(N+1,1);
dK2 = zeros(N+1,1);
K1 = zeros(N+1,1);
K2 = zeros(N+1,1);
dP = zeros(N+1,1);
P = zeros(N+1,1);
y = zeros(N+1,2);

y(1,1) = 1;
y(1,2) = 0.5;


for n=1:N

   k1 = h*f(y(n),t,K1,K2,dK1,dK2,P,dP);
   k2 = h*f(y(n)+1/2*k1,t,K1,K2,dK1,dK2,P,dP);
   k3 = h*f(y(n)+1/2*k2,t,K1,K2,dK1,dK2,P,dP);
   k4 = h*f(y(n)+k3,t,K1,K2,dK1,dK2,P,dP);

end

plot(t,y1,t,y2)

function dy = f(t,K1,K2,dK1,dK2,dP,P,y)

dy = zeros(1,2);

dy(1) = (dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^4 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^3 + 6.*dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^2 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2) + dP.*K1.^3.*K2.*P.^2.*y(1).^5 + dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^4 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^3 + 6.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^2 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2) + dP.*K1.^3.*K2.*P.^2.*y(1).^4 - 2.*dK2.*K1.^3.*P.^3.*y(1).^5.*y(2).^4 + 8.*dK2.*K1.^3.*P.^3.*y(1).^5.*y(2).^3 - 12.*dK2.*K1.^3.*P.^3.*y(1).^5.*y(2).^2 + 8.*dK2.*K1.^3.*P.^3.*y(1).^5.*y(2) - 2.*dK2.*K1.^3.*P.^3.*y(1).^5 - 2.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^4 + 8.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^3 - 12.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^2 + 8.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2) - 2.*dK2.*K1.^3.*P.^3.*y(1).^4 + 4.*dP.*K1.*K2.^3.*y(1).^5.*y(2).^2 + 8.*dP.*K1.*K2.^3.*y(1).^5.*y(2) + 4.*dP.*K1.*K2.^3.*y(1).^5 - 4.*dP.*K1.*K2.^3.*y(1).^4.*y(2).^2 + 8.*dP.*K1.*K2.^3.*y(1).^4.*y(2) + 12.*dP.*K1.*K2.^3.*y(1).^4 - 4.*dP.*K1.*K2.^3.*y(1).^3.*y(2).^2 - 24.*dP.*K1.*K2.^3.*y(1).^3.*y(2) - 4.*dP.*K1.*K2.^3.*y(1).^3 + 4.*dP.*K1.*K2.^3.*y(1).^2.*y(2).^2 - 8.*dP.*K1.*K2.^3.*y(1).^2.*y(2) - 28.*dP.*K1.*K2.^3.*y(1).^2 + 16.*dP.*K1.*K2.^3.*y(1).*y(2) + 16.*dP.*K1.*K2.^3 - 4.*dP.*K1.*K2.*P.*y(1).^5.*y(2).^3 + 4.*dP.*K1.*K2.*P.*y(1).^5.*y(2).^2 + 4.*dP.*K1.*K2.*P.*y(1).^5.*y(2) - 4.*dP.*K1.*K2.*P.*y(1).^5 + 8.*dP.*K1.*K2.*P.*y(1).^4.*y(2).^3 - 16.*dP.*K1.*K2.*P.*y(1).^4.*y(2).^2 + 8.*dP.*K1.*K2.*P.*y(1).^4.*y(2) - 4.*dP.*K1.*K2.*P.*y(1).^3.*y(2).^3 + 20.*dP.*K1.*K2.*P.*y(1).^3.*y(2).^2 - 28.*dP.*K1.*K2.*P.*y(1).^3.*y(2) + 12.*dP.*K1.*K2.*P.*y(1).^3 - 8.*dP.*K1.*K2.*P.*y(1).^2.*y(2).^2 + 16.*dP.*K1.*K2.*P.*y(1).^2.*y(2) - 8.*dP.*K1.*K2.*P.*y(1).^2 + 8.*dK2.*K1.*P.^2.*y(1).^5.*y(2).^3 + 8.*dK2.*K1.*P.^2.*y(1).^5.*y(2).^2 - 8.*dK2.*K1.*P.^2.*y(1).^5.*y(2) - 8.*dK2.*K1.*P.^2.*y(1).^5 - 16.*dK2.*K1.*P.^2.*y(1).^4.*y(2).^3 + 16.*dK2.*K1.*P.^2.*y(1).^4.*y(2).^2 + 16.*dK2.*K1.*P.^2.*y(1).^4.*y(2) - 16.*dK2.*K1.*P.^2.*y(1).^4 + 8.*dK2.*K1.*P.^2.*y(1).^3.*y(2).^3 - 56.*dK2.*K1.*P.^2.*y(1).^3.*y(2).^2 + 24.*dK2.*K1.*P.^2.*y(1).^3.*y(2) + 24.*dK2.*K1.*P.^2.*y(1).^3 + 32.*dK2.*K1.*P.^2.*y(1).^2.*y(2).^2 - 64.*dK2.*K1.*P.^2.*y(1).^2.*y(2) + 32.*dK2.*K1.*P.^2.*y(1).^2 + 32.*dK2.*K1.*P.^2.*y(1).*y(2) - 32.*dK2.*K1.*P.^2.*y(1) + 8.*dK1.*K2.^3.*P.*y(1).^5.*y(2).^2 + 16.*dK1.*K2.^3.*P.*y(1).^5.*y(2) + 8.*dK1.*K2.^3.*P.*y(1).^5 - 8.*dK1.*K2.^3.*P.*y(1).^4.*y(2).^2 + 16.*dK1.*K2.^3.*P.*y(1).^4.*y(2) + 24.*dK1.*K2.^3.*P.*y(1).^4 - 8.*dK1.*K2.^3.*P.*y(1).^3.*y(2).^2 - 48.*dK1.*K2.^3.*P.*y(1).^3.*y(2) - 8.*dK1.*K2.^3.*P.*y(1).^3 + 8.*dK1.*K2.^3.*P.*y(1).^2.*y(2).^2 - 16.*dK1.*K2.^3.*P.*y(1).^2.*y(2) - 56.*dK1.*K2.^3.*P.*y(1).^2 + 32.*dK1.*K2.^3.*P.*y(1).*y(2) + 32.*dK1.*K2.^3.*P + 16.*dK1.*K2.*P.^2.*y(1).^5.*y(2).^2 - 16.*dK1.*K2.*P.^2.*y(1).^5 - 16.*dK1.*K2.*P.^2.*y(1).^4.*y(2).^2 + 32.*dK1.*K2.*P.^2.*y(1).^4.*y(2) - 16.*dK1.*K2.*P.^2.*y(1).^4 - 16.*dK1.*K2.*P.^2.*y(1).^3.*y(2).^2 - 32.*dK1.*K2.*P.^2.*y(1).^3.*y(2) + 48.*dK1.*K2.*P.^2.*y(1).^3 + 16.*dK1.*K2.*P.^2.*y(1).^2.*y(2).^2 - 32.*dK1.*K2.*P.^2.*y(1).^2.*y(2) + 16.*dK1.*K2.*P.^2.*y(1).^2 + 32.*dK1.*K2.*P.^2.*y(1).*y(2) - 32.*dK1.*K2.*P.^2.*y(1))/(16.*K1.*K2.*P.*(y(1) - 1).*(y(1) + y(1).*y(2) + 2).*(K2.^2.*y(1) - 2.*P.*y(1) + 2.*K2.^2 + 2.*P.*y(1).*y(2) + K2.^2.*y(1).*y(2)));

dy(2) = -((y(2) - 1).*(dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^4 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^3 + 6.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2).^2 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^4.*y(2) + dP.*K1.^3.*K2.*P.^2.*y(1).^4 - 2.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^4 + 8.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^3 - 12.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2).^2 + 8.*dK2.*K1.^3.*P.^3.*y(1).^4.*y(2) - 2.*dK2.*K1.^3.*P.^3.*y(1).^4 + 4.*dP.*K1.*K2.^3.*y(1).^4.*y(2).^2 + 8.*dP.*K1.*K2.^3.*y(1).^4.*y(2) + 4.*dP.*K1.*K2.^3.*y(1).^4 - 8.*dP.*K1.*K2.^3.*y(1).^3.*y(2).^2 + 8.*dP.*K1.*K2.^3.*y(1).^3 + 4.*dP.*K1.*K2.^3.*y(1).^2.*y(2).^2 - 24.*dP.*K1.*K2.^3.*y(1).^2.*y(2) - 12.*dP.*K1.*K2.^3.*y(1).^2 + 16.*dP.*K1.*K2.^3.*y(1).*y(2) - 16.*dP.*K1.*K2.^3.*y(1) + 16.*dP.*K1.*K2.^3 + 4.*dP.*K1.*K2.*P.*y(1).^4.*y(2).^3 + 12.*dP.*K1.*K2.*P.*y(1).^4.*y(2).^2 - 4.*dP.*K1.*K2.*P.*y(1).^4.*y(2) - 12.*dP.*K1.*K2.*P.*y(1).^4 - 4.*dP.*K1.*K2.*P.*y(1).^3.*y(2).^3 - 4.*dP.*K1.*K2.*P.*y(1).^3.*y(2).^2 + 20.*dP.*K1.*K2.*P.*y(1).^3.*y(2) - 12.*dP.*K1.*K2.*P.*y(1).^3 - 8.*dP.*K1.*K2.*P.*y(1).^2.*y(2).^2 - 16.*dP.*K1.*K2.*P.*y(1).^2.*y(2) + 24.*dP.*K1.*K2.*P.*y(1).^2 - 8.*dK2.*K1.*P.^2.*y(1).^4.*y(2).^3 - 8.*dK2.*K1.*P.^2.*y(1).^4.*y(2).^2 + 8.*dK2.*K1.*P.^2.*y(1).^4.*y(2) + 8.*dK2.*K1.*P.^2.*y(1).^4 + 8.*dK2.*K1.*P.^2.*y(1).^3.*y(2).^3 - 24.*dK2.*K1.*P.^2.*y(1).^3.*y(2).^2 - 8.*dK2.*K1.*P.^2.*y(1).^3.*y(2) + 24.*dK2.*K1.*P.^2.*y(1).^3 + 32.*dK2.*K1.*P.^2.*y(1).^2.*y(2).^2 - 32.*dK2.*K1.*P.^2.*y(1).^2.*y(2) + 32.*dK2.*K1.*P.^2.*y(1).*y(2) - 32.*dK2.*K1.*P.^2.*y(1) + 8.*dK1.*K2.^3.*P.*y(1).^4.*y(2).^2 + 16.*dK1.*K2.^3.*P.*y(1).^4.*y(2) + 8.*dK1.*K2.^3.*P.*y(1).^4 - 16.*dK1.*K2.^3.*P.*y(1).^3.*y(2).^2 + 16.*dK1.*K2.^3.*P.*y(1).^3 + 8.*dK1.*K2.^3.*P.*y(1).^2.*y(2).^2 - 48.*dK1.*K2.^3.*P.*y(1).^2.*y(2) - 24.*dK1.*K2.^3.*P.*y(1).^2 + 32.*dK1.*K2.^3.*P.*y(1).*y(2) - 32.*dK1.*K2.^3.*P.*y(1) + 32.*dK1.*K2.^3.*P + 16.*dK1.*K2.*P.^2.*y(1).^4.*y(2).^2 - 16.*dK1.*K2.*P.^2.*y(1).^4 - 32.*dK1.*K2.*P.^2.*y(1).^3.*y(2).^2 + 32.*dK1.*K2.*P.^2.*y(1).^3.*y(2) + 16.*dK1.*K2.*P.^2.*y(1).^2.*y(2).^2 - 64.*dK1.*K2.*P.^2.*y(1).^2.*y(2) + 48.*dK1.*K2.*P.^2.*y(1).^2 + 32.*dK1.*K2.*P.^2.*y(1).*y(2) - 32.*dK1.*K2.*P.^2.*y(1)))/(16.*K1.*K2.*P.*y(1).*(y(1) - 1).*(y(1) + y(1).*y(2) + 2).*(K2.^2.*y(1) - 2.*P.*y(1) + 2.*K2.^2 + 2.*P.*y(1).*y(2) + K2.^2.*y(1).*y(2)));

dy = [dy(1); dy(2)];

end 

[Jan]

@Jason Louison: A general hint: Whenever you mention an error in the forum, post a copy of the complete message. Then the readers do not have to guess, which line is failing.

"Is there a way to do this in one step or command?" - seriously? The Replace dialog in the editor should do this sufficiently already. Even for really large functions with tenthousands of equations this takes fractions of a second only.

Using a fixed step Runge Kutta solver of order 4 solve a chemical reaction is a very bad idea. If the equation is stiff, the final value is dominated by the discretization error and more or less random. To reduce the noise, you have to run it with a tiny step size, such that the runtime will explode and the accumulated rounding errors will be the main part of the final value, such that you should not call it "result" anymore.

If the equation is stiff, you need a specific solver like ODE15s. Otherwise a solver with step-size control is the way to go, e.g. ODE45.

Solving an equation for a set of different parameters should not be done by vectorizing the function to be integrated, but by starting the integration in a loop, in which the parameters are varied. If multiple trajectories are calculated in the same integration, the component with the highest discretization and rounding errors determine the step-size. This means, that all other components are calculated too frequently and this reduces the accuracy due to the accumulation of rounding errors. The total run-time will increase because the system is integrated with a too small step-size in average.

If you have a really good reason to integrate a system with parameters provided as vectors, the variables to be integrated are not scalars anymore:

dy(1) = (dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^4 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^3 + ...

dy(1) is a scalar, but the right hand side is a vector. Of course this cannot work. You have to blow up dy to a 2*(N+1) vector instead.

I still do not see the reason to vectorize the equations. This does not offer an advantage for integrating the system. The shown Runge-Kutta method is not sufficient to calculate the trajectories reliably in addition. This is not state of the art even in lessons for beginners, but an expensive random number generator with a poor entropy.

Use on of the built-in solvers with a step-size control and check the stiffness of the system (a cheap check: Try, if ODE45 or ODE15s is much faster).

[Jason Louison]

I figured out a way to vectorize the function via the strrep function. That issue has been resolved. However, the reason as to why I asked how to do so in the first place was because the program kept giving me an “incorrect dimension error”. It was telling me that the operators I was using for the differential equation were incorrect, along the lines of

“Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform elementwise multiplication, use '.*'.”

No matter how hard I tried to fix the issue, it persisted. I was trying to avoid having to vectorize the differential equations in the first place.

Would you be willing to show some examples of code that follows exactly what you are suggesting in terms of using a loop to integrate? I am trying to teach myself to use matlab and/or simulink for my eventual internal combustion engine simulation program, of which this code is to function as a proof of concept for the modeling of combustion and dissociation within the cylinder, a very important part of the simulation.

P.S: No need for the harsh comments. I have taken the time here to become vulnerable and ask for help, which I usually never do. I am here to learn and better myself for both my academics and future, and we also all have to start somewhere. Nevertheless, Thank you for your help so far.

[Jan]

@Jason Louison: Again: "the program kept giving me an “incorrect dimension error”." - If you mention an error, post a copy of the complete message and the relevant lines of code.

My comments are not harsh. I try to show you, how to solve the problem and how to avoid the typical errors made by beginners. Numerical maths is a demanding field of science. Feel free to ask questions about Matlab, because this is the purpose of this forum.

What is the actual problem you want to solve? You have a 2 DoF system of equations. Why do you want to solve the system "as an array"? What do you expect from vectorizing the function to be integrated?

One of your posted codes contain:

dK = zeros(1,N+1);

dP = zeros(1,N+1);

P = zeros(1,N+1);

K = zeros(1,N+1);

%This allows us to vary the above parameters as functions of time. The

%parameters are assigned constant values below

dK = 0.34;

dP = 1;

P = 1;

K = 1;

After creating dK etc. as vectors, you overwrite all of them by scalars. What is the underlying intention?

An example for a proper integration:

dK1 = 0;
dK2 = 0;
K1 = 0.5;
K2 = 0.7; 
dP = 0;
P = 1;
[T, Y] = ode45(@(t,y) f(t,y, K1,K2,dK1,dK2,dP,P), [0, 60], [1, 0.5]);
plot(T, Y)

function dy = f(t,y, K1,K2,dK1,dK2,dP,P)
y1_5 = y(1)^5;
y1_4 = y(1)^4;
y1_3 = y(1)^3;
y1_2 = y(1)^2;

y2_4 = y(2)^4;
y2_3 = y(2)^3;
y2_2 = y(2)^2;

K2_3 = K2^3;
K2_2 = K2^2;
P_2  = P^2;

dy = [(dP*K1^3*K2*P_2*y1_5*y2_4 - 4*dP*K1^3*K2*P_2*y1_5*y2_3 + ...
       6*dP*K1^3*K2*P_2*y1_5*y2_2 - 4*dP*K1^3*K2*P_2*y1_5*y(2) + ...
       dP*K1^3*K2*P_2*y1_5 + dP*K1^3*K2*P_2*y1_4*y2_4 - ...
       4*dP*K1^3*K2*P_2*y1_4*y2_3 + 6*dP*K1^3*K2*P_2*y1_4*y2_2 - ...
       4*dP*K1^3*K2*P_2*y1_4*y(2) + dP*K1^3*K2*P_2*y1_4 - ...
       2*dK2*K1^3*P^3*y1_5*y2_4 + 8*dK2*K1^3*P^3*y1_5*y2_3 - ...
       12*dK2*K1^3*P^3*y1_5*y2_2 + 8*dK2*K1^3*P^3*y1_5*y(2) - ...
       2*dK2*K1^3*P^3*y1_5 - 2*dK2*K1^3*P^3*y1_4*y2_4 + ...
       8*dK2*K1^3*P^3*y1_4*y2_3 - 12*dK2*K1^3*P^3*y1_4*y2_2 + ...
       8*dK2*K1^3*P^3*y1_4*y(2) - 2*dK2*K1^3*P^3*y1_4 + ...
       4*dP*K1*K2_3*y1_5*y2_2 + 8*dP*K1*K2_3*y1_5*y(2) + ...
       4*dP*K1*K2_3*y1_5 - 4*dP*K1*K2_3*y1_4*y2_2 + ...
       8*dP*K1*K2_3*y1_4*y(2) + 12*dP*K1*K2_3*y1_4 - ...
       4*dP*K1*K2_3*y1_3*y2_2 - 24*dP*K1*K2_3*y1_3*y(2) - ...
       4*dP*K1*K2_3*y1_3 + 4*dP*K1*K2_3*y1_2*y2_2 - ...
       8*dP*K1*K2_3*y1_2*y(2) - 28*dP*K1*K2_3*y1_2 + ...
       16*dP*K1*K2_3*y(1)*y(2) + 16*dP*K1*K2_3 - ...
       4*dP*K1*K2*P*y1_5*y2_3 + 4*dP*K1*K2*P*y1_5*y2_2 + ...
       4*dP*K1*K2*P*y1_5*y(2) - 4*dP*K1*K2*P*y1_5 + ...
       8*dP*K1*K2*P*y1_4*y2_3 - 16*dP*K1*K2*P*y1_4*y2_2 + ...
       8*dP*K1*K2*P*y1_4*y(2) - 4*dP*K1*K2*P*y1_3*y2_3 + ...
       20*dP*K1*K2*P*y1_3*y2_2 - 28*dP*K1*K2*P*y1_3*y(2) + ...
       12*dP*K1*K2*P*y1_3 - 8*dP*K1*K2*P*y1_2*y2_2 + ...
       16*dP*K1*K2*P*y1_2*y(2) - 8*dP*K1*K2*P*y1_2 + ...
       8*dK2*K1*P_2*y1_5*y2_3 + 8*dK2*K1*P_2*y1_5*y2_2 - ...
       8*dK2*K1*P_2*y1_5*y(2) - 8*dK2*K1*P_2*y1_5 - ...
       16*dK2*K1*P_2*y1_4*y2_3 + 16*dK2*K1*P_2*y1_4*y2_2 + ...
       16*dK2*K1*P_2*y1_4*y(2) - 16*dK2*K1*P_2*y1_4 + ...
       8*dK2*K1*P_2*y1_3*y2_3 - 56*dK2*K1*P_2*y1_3*y2_2 + ...
       24*dK2*K1*P_2*y1_3*y(2) + 24*dK2*K1*P_2*y1_3 + ...
       32*dK2*K1*P_2*y1_2*y2_2 - 64*dK2*K1*P_2*y1_2*y(2) + ...
       32*dK2*K1*P_2*y1_2 + 32*dK2*K1*P_2*y(1)*y(2) - ...
       32*dK2*K1*P_2*y(1) + 8*dK1*K2_3*P*y1_5*y2_2 + ...
       16*dK1*K2_3*P*y1_5*y(2) + 8*dK1*K2_3*P*y1_5 - ...
       8*dK1*K2_3*P*y1_4*y2_2 + 16*dK1*K2_3*P*y1_4*y(2) + ...
       24*dK1*K2_3*P*y1_4 - 8*dK1*K2_3*P*y1_3*y2_2 - ...
       48*dK1*K2_3*P*y1_3*y(2) - 8*dK1*K2_3*P*y1_3 + ...
       8*dK1*K2_3*P*y1_2*y2_2 - 16*dK1*K2_3*P*y1_2*y(2) - ...
       56*dK1*K2_3*P*y1_2 + 32*dK1*K2_3*P*y(1)*y(2) + ...
       32*dK1*K2_3*P + 16*dK1*K2*P_2*y1_5*y2_2 - ...
       16*dK1*K2*P_2*y1_5 - 16*dK1*K2*P_2*y1_4*y2_2 + ...
       32*dK1*K2*P_2*y1_4*y(2) - 16*dK1*K2*P_2*y1_4 - ...
       16*dK1*K2*P_2*y1_3*y2_2 - 32*dK1*K2*P_2*y1_3*y(2) + ...
       48*dK1*K2*P_2*y1_3 + 16*dK1*K2*P_2*y1_2*y2_2 - ...
       32*dK1*K2*P_2*y1_2*y(2) + 16*dK1*K2*P_2*y1_2 + ...
       32*dK1*K2*P_2*y(1)*y(2) - 32*dK1*K2*P_2*y(1)) / ...
       (16*K1*K2*P*(y(1) - 1)*(y(1) + y(1)*y(2) + 2)*(K2_2*y(1) - ...
       2*P*y(1) + 2*K2_2 + 2*P*y(1)*y(2) + K2_2*y(1)*y(2))); ...
     -((y(2) - 1)*(dP*K1^3*K2*P_2*y1_4*y2_4 - 4*dP*K1^3*K2*P_2*y1_4*y2_3 + ...
	 6*dP*K1^3*K2*P_2*y1_4*y2_2 - 4*dP*K1^3*K2*P_2*y1_4*y(2) + ...
	 dP*K1^3*K2*P_2*y1_4 - 2*dK2*K1^3*P^3*y1_4*y2_4 + ...
	 8*dK2*K1^3*P^3*y1_4*y2_3 - 12*dK2*K1^3*P^3*y1_4*y2_2 + ...
	 8*dK2*K1^3*P^3*y1_4*y(2) - 2*dK2*K1^3*P^3*y1_4 + ...
	 4*dP*K1*K2_3*y1_4*y2_2 + 8*dP*K1*K2_3*y1_4*y(2) + ...
	 4*dP*K1*K2_3*y1_4 - 8*dP*K1*K2_3*y1_3*y2_2 + ...
	 8*dP*K1*K2_3*y1_3 + 4*dP*K1*K2_3*y1_2*y2_2 - ...
	 24*dP*K1*K2_3*y1_2*y(2) - 12*dP*K1*K2_3*y1_2 + ...
	 16*dP*K1*K2_3*y(1)*y(2) - 16*dP*K1*K2_3*y(1) + ...
	 16*dP*K1*K2_3 + 4*dP*K1*K2*P*y1_4*y2_3 + ...
	 12*dP*K1*K2*P*y1_4*y2_2 - 4*dP*K1*K2*P*y1_4*y(2) - ...
	 12*dP*K1*K2*P*y1_4 - 4*dP*K1*K2*P*y1_3*y2_3 - ...
	 4*dP*K1*K2*P*y1_3*y2_2 + 20*dP*K1*K2*P*y1_3*y(2) - ...
	 12*dP*K1*K2*P*y1_3 - 8*dP*K1*K2*P*y1_2*y2_2 - ...
	 16*dP*K1*K2*P*y1_2*y(2) + 24*dP*K1*K2*P*y1_2 - ...
	 8*dK2*K1*P_2*y1_4*y2_3 - 8*dK2*K1*P_2*y1_4*y2_2 + ...
	 8*dK2*K1*P_2*y1_4*y(2) + 8*dK2*K1*P_2*y1_4 + ...
	 8*dK2*K1*P_2*y1_3*y2_3 - 24*dK2*K1*P_2*y1_3*y2_2 - ...
	 8*dK2*K1*P_2*y1_3*y(2) + 24*dK2*K1*P_2*y1_3 + ...
	 32*dK2*K1*P_2*y1_2*y2_2 - 32*dK2*K1*P_2*y1_2*y(2) + ...
	 32*dK2*K1*P_2*y(1)*y(2) - 32*dK2*K1*P_2*y(1) + ...
	 8*dK1*K2_3*P*y1_4*y2_2 + 16*dK1*K2_3*P*y1_4*y(2) + ...
	 8*dK1*K2_3*P*y1_4 - 16*dK1*K2_3*P*y1_3*y2_2 + ...
	 16*dK1*K2_3*P*y1_3 + 8*dK1*K2_3*P*y1_2*y2_2 - ...
	 48*dK1*K2_3*P*y1_2*y(2) - 24*dK1*K2_3*P*y1_2 + ...
	 32*dK1*K2_3*P*y(1)*y(2) - 32*dK1*K2_3*P*y(1) + ...
	 32*dK1*K2_3*P + 16*dK1*K2*P_2*y1_4*y2_2 - ...
	 16*dK1*K2*P_2*y1_4 - 32*dK1*K2*P_2*y1_3*y2_2 + ...
	 32*dK1*K2*P_2*y1_3*y(2) + 16*dK1*K2*P_2*y1_2*y2_2 - ...
	 64*dK1*K2*P_2*y1_2*y(2) + 48*dK1*K2*P_2*y1_2 + ...
	 32*dK1*K2*P_2*y(1)*y(2) - 32*dK1*K2*P_2*y(1))) / ...
	 (16*K1*K2*P*y(1)*(y(1) - 1)*(y(1) + y(1)*y(2) + 2)*(K2_2*y(1) - ...
	 2*P*y(1) + 2*K2_2 + 2*P*y(1)*y(2) + K2_2*y(1)*y(2)))];
end

The output is empty, because the calculated trajectory contains NaN values only. I assume it is divided by zero.

So a vectorization is not meaningful at all. It is not clear yet, what you want to achieve and the posted functions does not have a useful result.

[Jason Louison]

The reason I overwrite the values as scalars is because I do not know how to make each value in the matrix assigned to each array a constant value. For example, I want each value of “dK” to be 1 for all time t. As I do not yet know how these values will actually vary with time in the real simulation, I would like the option to switch from constant for all time to to being able to write each parameter as a function of time t. Both in array format.

Answer 1

[Torsten]

See how you can use the below code for your purpose.

The Runge-Kutta code is made to deal with systems of differential equations of arbitrary size (thus with vector inputs).

tstart = 0.0;
tend = 1.0;
nt = 100;
h = (tend - tstart)/nt;
T = (tstart:h:tend).';
Y0 = [1 -1]; 
B = 4;
Y = runge_kutta_RK4(@(t,y)f(t,y,B),T,Y0);
plot(T,Y)                          

function Y = runge_kutta_RK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    t = T(i-1);
    y = Y(i-1,:);
    h  = T(i) - T(i-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end 

function dy = f(t,y,B)
  dy = zeros(1,2);
  dy(1) = y(2);
  dy(2) = -exp(-B*t)-y(1)+5*exp(-2*t)-2*exp(-(B+2)*t)+exp(-B*t)+t;
end

[Jason Louison]

I will try implementing this code when I get a chance today. Will update when I do so.

[Jason Louison]

Tortsen,

After hours of trial and error and research, I have finally got the program to work. It does exactly what I want it to do. I cannot thank you and Jan enough. The latest version of the code is shown below. Attached is the output for the given input parameters and initial values.

I do believe the program can be optimized even more, especially with the use of the ode45 solver, but I wanted to see a runge kutta 4th order algorithm be put to use. Other areas of improvement could be a more organized layout as well as a less confusing variable input system (The functions MUST have their own variables, and the values I assign within the function don't get used. I find it strange). Another worry of mine is that the initial conditions cannot be zero for both y(1) and y(2). One has to be greater than zero, and calculating this right before the reaction begins may prove difficult. The main flaw of the program that was keeping me from being able to solve the system was that I needed a for loop to tell the differential equations what the value of each parameter was each time step. Nevertheless, I am extremely happy and grateful, as this code serves as a basis for solving more complex reactions. Here is a result figure I obtained from a different run, now showing the amount of each substance as well as dissociation and mole fraction:

% Chemical Equilibrium

% Solve the following chemical equation: 
% aCO2 -> n1[CO2]+n2[O2]+n3[CO]+n4[O]

% Goal: We want to be able to solve the chemical reaction above for all
% products assuming chemical equilibrium

% Assume a fraction of CO2 and O2 dissociates. We will call these parameters
% "y1" and "y2" 

% Write out the mass balance: 
% a -> a(1-y1)[CO2]+(1/2)ay1(1-y2)[O2]+ay1[CO]+ay1[CO]+ay1y2[O]

% From this, we can derive the total number of products:

% np = a[1-y1+(1/2)y1(1-y2)+y1+y1y2]
% np = (1/2)a[2+y1(1+y2)]

% From this, we define each mole fraction

% x1 = 2(1-y1)/(2+y1(1+y2))
% x2 = y1(1-y2)/(2+y1(1+y2))
% x3 = 2y1/(2+y1(1+y2))
% x4 = 2y1y2/(2+y1(1+y2))

% From here, we can plug these values into the equilibrium constant
% equation and solve for the time derivate of y. To solve for the
% differential equation, initialization parameters xij are used.

%clearvars

% Initialization Parameters: These mole fraction and dissociation constant
% equations are being used to solve for the time derivatives of y1 and y2

clearvars

syms K1 dK1 K2 dK2 y1 dy1 y2 dy2 P dP 

x1j = 2*(1-y1)/(2+y1*(1+y2));

x2j = y1*(1-y2)/(2+y1*(1+y2));

x3j = 2*y1/(2+y1*(1+y2));

x4j = (2*y1*y2)/(2+y1*(1+y2));

% Differential Initialization Parameters

dx1j = diff(x1j,y1)*dy1+diff(x1j,y2)*dy2;

dx2j = diff(x2j,y1)*dy1+diff(x2j,y2)*dy2;

dx3j = diff(x3j,y1)*dy1+diff(x3j,y2)*dy2;

dx4j = diff(x4j,y1)*dy1+diff(x4j,y2)*dy2;

% Dissociation Constant Initialization Parameters

x2jk = (x4j^2*P)/(K2^2);

x3jk = (x1j^2)/(x2j^2*K1^2*P);

%solving the combined equations for dy1 and dy2

eqn = [(diff(x2j,y1)*dy1 + diff(x2j,y2)*dy2 == diff(x2jk,y1)*dy1 + diff(x2jk,y2)*dy2 + diff(x2jk,P)*dP + diff(x2jk,K1)*dK1 + diff(x2jk,K2)*dK2)
    (diff(x3j,y1)*dy1 + diff(x3j,y2)*dy2 == diff(x3jk,y1)*dy1 + diff(x3jk,y2)*dy2 + diff(x3jk,P)*dP + diff(x3jk,K1)*dK1 + diff(x3jk,K2)*dK2)];

sol = solve(eqn,[dy1 dy2]);

dy1i = simplify(sol.dy1,"Steps",50);

dy2i = simplify(sol.dy2,"Steps",50);


%Copy and Paste the equations from the variables dy1i and dy2i into the
%differential equation functions at the bottom of the code. Be sure to
%replace all operators with array compatible ones, all parametes with a t
%at the end, and y1 and y2 to y(1) and y(2)


%Now that we have equations for dy1 and dy2, we can proceed with the
%numerical analysis of the system



clearvars

tstart = 0.0;
tend = 120.0;
nt = 1000000;
h = (tend - tstart)/nt;
T = (tstart:h:tend).';

%Values of Parameters. These are stored in the matrix elements for each. It
%can be a function of time or a constant value. 

%Global Variables

dK1 = sin(3.14*T);
dK2 = 2*cos(3.14*T);
K1 = 60;
K2 = 25;
dP = 1;
P = 1;

%Initial Conditions for Dissociation fraction y1 and y2

Y0 = [0.01 0.00];  


%Array variables defining global varables at time t as an array

dK1f = dK1 + zeros(nt+1,1);
dK2f = dK2 + zeros(nt+1,1);
K1f = K1 + zeros(nt+1,1);
K2f = K2 + zeros(nt+1,1);
dPf = dP + zeros(nt+1,1);
Pf = P + zeros(nt+1,1);
yf = zeros(nt+1,2);

N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;

  for i = 2:N
  dK1t = dK1f(i-1);
  dK2t = dK2f(i-1);
  K1t = K1f(i-1);
  K2t = K2f(i-1);
  dPt = dPf(i-1);
  Pt = Pf(i-1);
  end

Y = runge_kutta_RK4(@(t,y)f(t,K1t,K2t,dK1t,dK2t,dPt,Pt,y),T,Y0);

x1 = 2.*(1-Y(:,1))./(2+Y(:,1).*(1+Y(:,2)));

x2 = (Y(:,1).*(1-Y(:,2)))./(2+Y(:,1).*(1+Y(:,2)));

x3 = 2.*Y(:,1)./(2+Y(:,1).*(1+Y(:,2)));

x4 = (2.*Y(:,1).*Y(:,2))./(2+Y(:,1).*(1+Y(:,2)));


subplot(1,2,1);
plot(T,Y(:,1),T,Y(:,2))
title('Dissociation Fraction of CO_2 and O_2 vs. Time')
xlabel('Time (s)')
ylabel('Dissociation Fraction (ni/nt)')
legend({'CO_2','O_2'})
subplot(1,2,2)
plot(T,x1,T,x2,T,x3,T,x4)
title('Mole Fractions of CO_2, O_2, CO, and O vs. Time')
xlabel('Time (s)')
ylabel('Mole Fraction (ni/nt)')
legend({'CO_2','O_2','CO','O'})

% Differential Equation Functions 

function dy = f(~,dK1t,dK2t,K1t,K2t,dPt,Pt,y)

dy = zeros(1,2);

dy(1) = -(2.*(- dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^5 + 3.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^4 - 3.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^3 + ...
    dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^2 + 2.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^5 - 6.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^4 + 6.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^3 - ...
    2.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^5.*y(2).^3 - dPt.*K1t.*K2t.^3.*y(1).^5.*y(2).^2 + dPt.*K1t.*K2t.^3.*y(1).^5.*y(2) + dPt.*K1t.*K2t.^3.*y(1).^5 + ...
    dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^3 - 3.*dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^4.*y(2) + 3.*dPt.*K1t.*K2t.^3.*y(1).^4 + dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^3 + ...
    5.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^2 - 5.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2) - dPt.*K1t.*K2t.^3.*y(1).^3 - dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^3 + 3.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^2 + ...
    5.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2) - 7.*dPt.*K1t.*K2t.^3.*y(1).^2 - 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2).^2 + 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2) - 4.*dPt.*K1t.*K2t.^3.*y(2) + ...
    4.*dPt.*K1t.*K2t.^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^4 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^2 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2) - ...
    8.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^4 + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2) + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^4 - 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^3 - ...
    12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^2 - 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2) + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^3 + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^2 - ...
    16.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2) + 16.*dPt.*K1t.*K2t.*Pt.*y(1).*y(2) - 8.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^3 - ...
    8.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^2 + 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^2 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^4 + ...
    48.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^3 + 24.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^2 - 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^3 + 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^2 - ...
    32.*dK2t.*K1t.*Pt.^2.*y(1).*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2).^3 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2).^2 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2) + ...
    2.*dK1t.*K2t.^3.*Pt.*y(1).^5 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^3 - 6.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2) + 6.*dK1t.*K2t.^3.*Pt.*y(1).^4 + ...
    2.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^3 + 10.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^2 - 10.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2) - 2.*dK1t.*K2t.^3.*Pt.*y(1).^3 - ...
    2.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^3 + 6.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^2 + 10.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2) - 14.*dK1t.*K2t.^3.*Pt.*y(1).^2 - ...
    8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2).^2 + 8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2) - 8.*dK1t.*K2t.^3.*Pt.*y(2) + 8.*dK1t.*K2t.^3.*Pt - 8.*dK1t.*K2t.*Pt.^2.*y(1).^5.*y(2).^3 + ...
    8.*dK1t.*K2t.*Pt.^2.*y(1).^5.*y(2) - 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2) + 24.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2).^3 - ...
    24.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2) - 16.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^3 + 48.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^2 - 32.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2) - ...
    32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2).^2 + 32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2)))/(K1t.*K2t.*Pt.*(- K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^4 + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^3 - ...
    2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2) + K1t.^2.*K2t.^2.*Pt.*y(1).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^3 + 6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^2 - ...
    6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2) + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2 - 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^4 + 24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^3 - ...
    24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^2 + 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2) + 8.*K2t.^2.*y(1).^3.*y(2).^3 + 8.*K2t.^2.*y(1).^3.*y(2).^2 - 8.*K2t.^2.*y(1).^3.*y(2) - ...
    8.*K2t.^2.*y(1).^3 - 8.*K2t.^2.*y(1).^2.*y(2).^3 + 24.*K2t.^2.*y(1).^2.*y(2).^2 + 8.*K2t.^2.*y(1).^2.*y(2) - 24.*K2t.^2.*y(1).^2 - ...
    32.*K2t.^2.*y(1).*y(2).^2 + 32.*K2t.^2.*y(1).*y(2) - 32.*K2t.^2.*y(2) + 32.*K2t.^2 + 48.*Pt.*y(1).^3.*y(2).^3 + 32.*Pt.*y(1).^3.*y(2).^2 - 16.*Pt.*y(1).^3.*y(2) - ...
    48.*Pt.*y(1).^2.*y(2).^3 + 32.*Pt.*y(1).^2.*y(2).^2 - 48.*Pt.*y(1).^2.*y(2) - 64.*Pt.*y(1).*y(2).^2 + 64.*Pt.*y(2)));

dy(2) = -(2.*(y(2) - 1).*(- 2.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^4 + 4.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^3 - 2.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 + ...
    4.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^4 - 8.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^3 + 4.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^2 + dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^3 + ...
    dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^4.*y(2) - dPt.*K1t.*K2t.^3.*y(1).^4 - 2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^3 + 2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^2 + ...
    2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2) - 2.*dPt.*K1t.*K2t.^3.*y(1).^3 + dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^3 - 7.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^2 + ...
    3.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2) + 3.*dPt.*K1t.*K2t.^3.*y(1).^2 + 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2).^2 - 8.*dPt.*K1t.*K2t.^3.*y(1).*y(2) + 4.*dPt.*K1t.*K2t.^3.*y(1) + ...
    4.*dPt.*K1t.*K2t.^3.*y(2) - 4.*dPt.*K1t.*K2t.^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^4 + 16.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^3 + 12.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^2 - ...
    4.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^4 - 8.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^3 + 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^2 - 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^3 - ...
    24.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^2 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^3 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^2 + ...
    8.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^3 - 24.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^2 + 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^3 + ...
    32.*dK2t.*K1t.*Pt.^2.*y(1).*y(2).^2 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^3 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2) - ...
    2.*dK1t.*K2t.^3.*Pt.*y(1).^4 - 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^3 + 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^2 + 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2) - ...
    4.*dK1t.*K2t.^3.*Pt.*y(1).^3 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^3 - 14.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^2 + 6.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2) + ...
    6.*dK1t.*K2t.^3.*Pt.*y(1).^2 + 8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2).^2 - 16.*dK1t.*K2t.^3.*Pt.*y(1).*y(2) + 8.*dK1t.*K2t.^3.*Pt.*y(1) + 8.*dK1t.*K2t.^3.*Pt.*y(2) - ...
    8.*dK1t.*K2t.^3.*Pt + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^3 + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 - 32.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2).^3 + ...
    16.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^3 - 48.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^2 + 32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2).^2))/(K1t.*K2t.*Pt.*y(1).*(- K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^4 + ...
    2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2) + K1t.^2.*K2t.^2.*Pt.*y(1).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^3 + ...
    6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^2 - 6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2) + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2 - 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^4 + ...
    24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^3 - 24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^2 + 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2) + 8.*K2t.^2.*y(1).^3.*y(2).^3 + ...
    8.*K2t.^2.*y(1).^3.*y(2).^2 - 8.*K2t.^2.*y(1).^3.*y(2) - 8.*K2t.^2.*y(1).^3 - 8.*K2t.^2.*y(1).^2.*y(2).^3 + 24.*K2t.^2.*y(1).^2.*y(2).^2 + ...
    8.*K2t.^2.*y(1).^2.*y(2) - 24.*K2t.^2.*y(1).^2 - 32.*K2t.^2.*y(1).*y(2).^2 + 32.*K2t.^2.*y(1).*y(2) - 32.*K2t.^2.*y(2) + 32.*K2t.^2 + 48.*Pt.*y(1).^3.*y(2).^3 + ...
    32.*Pt.*y(1).^3.*y(2).^2 - 16.*Pt.*y(1).^3.*y(2) - 48.*Pt.*y(1).^2.*y(2).^3 + 32.*Pt.*y(1).^2.*y(2).^2 - 48.*Pt.*y(1).^2.*y(2) - 64.*Pt.*y(1).*y(2).^2 + 64.*Pt.*y(2)));

dy = [dy(1); dy(2)];

end 

function Y = runge_kutta_RK4(f,T,Y0)
 
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;

dK1f = 1;
dK2f = 1;
K1f = 1;
K2f = 1;
dPf = 1;
Pf = 1;

nt = 10000000;

dK1 = dK1f + zeros(nt+1,1);
dK2 = dK2f + zeros(nt+1,1);
K1 = K1f + zeros(nt+1,1);
K2 = K2f + zeros(nt+1,1);
dP = dPf + zeros(nt+1,1);
P = Pf + zeros(nt+1,1);
y = zeros(nt+1,2);

  for i = 2:N
    
   dK1t = dK1(i-1);
   dK2t = dK2(i-1);
   K1t = K1(i-1);
   K2t = K2(i-1);
   dPt = dP(i-1);
   Pt = P(i-1);
   t = T(i-1);
   y = Y(i-1,:);
   h  = T(i) - T(i-1);
   k0 = h*f(t,y)';
   k1 = h*f(t,y+k0*0.5)';
   k2 = h*f(t,y+k1*0.5)';
   k3 = h*f(t,y+k2)';
   Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  
  end
end 

[Torsten]

I get different results.

I wonder why you modified the function "runge_kutta_RK4" and made it work incorrectly. Changes to "runge_kutta_RK4" are never necessary. Changes should to be made in the calling function and in the function where you supply the derivatives. Or do you modify MATLAB's inbuilt ODE45 when you want to solve a differential equation ?

tstart = 0.0;
tend = 120.0;
nt = 1000000;
h = (tend - tstart)/nt;
T = (tstart:h:tend).';
%Initial Conditions for Dissociation fraction y1 and y2
Y0 = [0.01 0.00];  
N = numel(T);
n = numel(Y0);
Y = runge_kutta_RK4(@f,T,Y0);
x1 = 2.*(1-Y(:,1))./(2+Y(:,1).*(1+Y(:,2)));
x2 = (Y(:,1).*(1-Y(:,2)))./(2+Y(:,1).*(1+Y(:,2)));
x3 = 2.*Y(:,1)./(2+Y(:,1).*(1+Y(:,2)));
x4 = (2.*Y(:,1).*Y(:,2))./(2+Y(:,1).*(1+Y(:,2)));
subplot(1,2,1);
plot(T,Y(:,1),T,Y(:,2))
title('Dissociation Fraction of CO_2 and O_2 vs. Time')
xlabel('Time (s)')
ylabel('Dissociation Fraction (ni/nt)')
legend({'CO_2','O_2'})
subplot(1,2,2)
plot(T,x1,T,x2,T,x3,T,x4)
title('Mole Fractions of CO_2, O_2, CO, and O vs. Time')
xlabel('Time (s)')
ylabel('Mole Fraction (ni/nt)')
legend({'CO_2','O_2','CO','O'})

% Differential Equation Functions 
function dy = f(t,y)
dK1t = sin(3.14*t);
dK2t = 2*cos(3.14*t);
K1t = 60;
K2t = 25;
dPt = 1;
Pt = 1;
dy = zeros(1,2);
dy(1) = -(2.*(- dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^5 + 3.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^4 - 3.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^3 + ...
    dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^2 + 2.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^5 - 6.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^4 + 6.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^3 - ...
    2.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^5.*y(2).^3 - dPt.*K1t.*K2t.^3.*y(1).^5.*y(2).^2 + dPt.*K1t.*K2t.^3.*y(1).^5.*y(2) + dPt.*K1t.*K2t.^3.*y(1).^5 + ...
    dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^3 - 3.*dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^4.*y(2) + 3.*dPt.*K1t.*K2t.^3.*y(1).^4 + dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^3 + ...
    5.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^2 - 5.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2) - dPt.*K1t.*K2t.^3.*y(1).^3 - dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^3 + 3.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^2 + ...
    5.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2) - 7.*dPt.*K1t.*K2t.^3.*y(1).^2 - 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2).^2 + 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2) - 4.*dPt.*K1t.*K2t.^3.*y(2) + ...
    4.*dPt.*K1t.*K2t.^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^4 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^2 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2) - ...
    8.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^4 + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2) + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^4 - 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^3 - ...
    12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^2 - 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2) + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^3 + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^2 - ...
    16.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2) + 16.*dPt.*K1t.*K2t.*Pt.*y(1).*y(2) - 8.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^3 - ...
    8.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^2 + 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^2 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^4 + ...
    48.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^3 + 24.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^2 - 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^3 + 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^2 - ...
    32.*dK2t.*K1t.*Pt.^2.*y(1).*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2).^3 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2).^2 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2) + ...
    2.*dK1t.*K2t.^3.*Pt.*y(1).^5 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^3 - 6.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2) + 6.*dK1t.*K2t.^3.*Pt.*y(1).^4 + ...
    2.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^3 + 10.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^2 - 10.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2) - 2.*dK1t.*K2t.^3.*Pt.*y(1).^3 - ...
    2.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^3 + 6.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^2 + 10.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2) - 14.*dK1t.*K2t.^3.*Pt.*y(1).^2 - ...
    8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2).^2 + 8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2) - 8.*dK1t.*K2t.^3.*Pt.*y(2) + 8.*dK1t.*K2t.^3.*Pt - 8.*dK1t.*K2t.*Pt.^2.*y(1).^5.*y(2).^3 + ...
    8.*dK1t.*K2t.*Pt.^2.*y(1).^5.*y(2) - 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2) + 24.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2).^3 - ...
    24.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2) - 16.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^3 + 48.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^2 - 32.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2) - ...
    32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2).^2 + 32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2)))/(K1t.*K2t.*Pt.*(- K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^4 + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^3 - ...
    2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2) + K1t.^2.*K2t.^2.*Pt.*y(1).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^3 + 6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^2 - ...
    6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2) + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2 - 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^4 + 24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^3 - ...
    24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^2 + 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2) + 8.*K2t.^2.*y(1).^3.*y(2).^3 + 8.*K2t.^2.*y(1).^3.*y(2).^2 - 8.*K2t.^2.*y(1).^3.*y(2) - ...
    8.*K2t.^2.*y(1).^3 - 8.*K2t.^2.*y(1).^2.*y(2).^3 + 24.*K2t.^2.*y(1).^2.*y(2).^2 + 8.*K2t.^2.*y(1).^2.*y(2) - 24.*K2t.^2.*y(1).^2 - ...
    32.*K2t.^2.*y(1).*y(2).^2 + 32.*K2t.^2.*y(1).*y(2) - 32.*K2t.^2.*y(2) + 32.*K2t.^2 + 48.*Pt.*y(1).^3.*y(2).^3 + 32.*Pt.*y(1).^3.*y(2).^2 - 16.*Pt.*y(1).^3.*y(2) - ...
    48.*Pt.*y(1).^2.*y(2).^3 + 32.*Pt.*y(1).^2.*y(2).^2 - 48.*Pt.*y(1).^2.*y(2) - 64.*Pt.*y(1).*y(2).^2 + 64.*Pt.*y(2)));
dy(2) = -(2.*(y(2) - 1).*(- 2.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^4 + 4.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^3 - 2.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 + ...
    4.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^4 - 8.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^3 + 4.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^2 + dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^3 + ...
    dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^4.*y(2) - dPt.*K1t.*K2t.^3.*y(1).^4 - 2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^3 + 2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^2 + ...
    2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2) - 2.*dPt.*K1t.*K2t.^3.*y(1).^3 + dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^3 - 7.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^2 + ...
    3.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2) + 3.*dPt.*K1t.*K2t.^3.*y(1).^2 + 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2).^2 - 8.*dPt.*K1t.*K2t.^3.*y(1).*y(2) + 4.*dPt.*K1t.*K2t.^3.*y(1) + ...
    4.*dPt.*K1t.*K2t.^3.*y(2) - 4.*dPt.*K1t.*K2t.^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^4 + 16.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^3 + 12.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^2 - ...
    4.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^4 - 8.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^3 + 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^2 - 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^3 - ...
    24.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^2 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^3 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^2 + ...
    8.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^3 - 24.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^2 + 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^3 + ...
    32.*dK2t.*K1t.*Pt.^2.*y(1).*y(2).^2 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^3 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2) - ...
    2.*dK1t.*K2t.^3.*Pt.*y(1).^4 - 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^3 + 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^2 + 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2) - ...
    4.*dK1t.*K2t.^3.*Pt.*y(1).^3 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^3 - 14.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^2 + 6.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2) + ...
    6.*dK1t.*K2t.^3.*Pt.*y(1).^2 + 8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2).^2 - 16.*dK1t.*K2t.^3.*Pt.*y(1).*y(2) + 8.*dK1t.*K2t.^3.*Pt.*y(1) + 8.*dK1t.*K2t.^3.*Pt.*y(2) - ...
    8.*dK1t.*K2t.^3.*Pt + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^3 + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 - 32.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2).^3 + ...
    16.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^3 - 48.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^2 + 32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2).^2))/(K1t.*K2t.*Pt.*y(1).*(- K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^4 + ...
    2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2) + K1t.^2.*K2t.^2.*Pt.*y(1).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^3 + ...
    6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^2 - 6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2) + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2 - 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^4 + ...
    24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^3 - 24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^2 + 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2) + 8.*K2t.^2.*y(1).^3.*y(2).^3 + ...
    8.*K2t.^2.*y(1).^3.*y(2).^2 - 8.*K2t.^2.*y(1).^3.*y(2) - 8.*K2t.^2.*y(1).^3 - 8.*K2t.^2.*y(1).^2.*y(2).^3 + 24.*K2t.^2.*y(1).^2.*y(2).^2 + ...
    8.*K2t.^2.*y(1).^2.*y(2) - 24.*K2t.^2.*y(1).^2 - 32.*K2t.^2.*y(1).*y(2).^2 + 32.*K2t.^2.*y(1).*y(2) - 32.*K2t.^2.*y(2) + 32.*K2t.^2 + 48.*Pt.*y(1).^3.*y(2).^3 + ...
    32.*Pt.*y(1).^3.*y(2).^2 - 16.*Pt.*y(1).^3.*y(2) - 48.*Pt.*y(1).^2.*y(2).^3 + 32.*Pt.*y(1).^2.*y(2).^2 - 48.*Pt.*y(1).^2.*y(2) - 64.*Pt.*y(1).*y(2).^2 + 64.*Pt.*y(2)));
end 
function Y = runge_kutta_RK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    t = T(i-1);
    y = Y(i-1,:);
    h  = T(i) - T(i-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end 

[Jason Louison]

Your differential equation was different, and my parameters had to be modified to change with time. This means that for each time step, a different value of dK1, dK2, dP…etc has to be used. I made a for loop to calculate or “find” the values of all the parameters at a certain time “t”. Yes, for some reason, your mole fraction is above 1. That is incorrect. I don’t know why your solution didn’t converge properly, but my code will not run properly unless it is exactly like the one I posted. I am a bit tried for time at the moment, but try running the exact code I had success with. Sometimes adding more time steps solves the issue of solution divergence, but it could also be your initial values and parameters.

[Torsten]

As said, already the changes made to the RK integration scheme made your code wrong:

k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);

is not the same as

k0 = h*f(t,y)';
k1 = h*f(t,y+k0*0.5)';
k2 = h*f(t,y+k1*0.5)';
k3 = h*f(t,y+k2)';
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);

And the time-dependent parameters are incorporated in f:

dK1t = sin(3.14*t);
dK2t = 2*cos(3.14*t);
K1t = 60;
K2t = 25;
dPt = 1;
Pt = 1;

[Jason Louison]

I had to do the for loop for the time dependent parameters because my code would not work otherwise. It returns an array size mismatch error. As for the runge kutta algorithm, the function itself, the differential equations, aren’t functions of time, but of y(1), and y(2). That is why it works. I have done the same with a spring mass damper simulation. My question to you is, would adding back the time parameter to the algorithm influence the solution behavior?

[Torsten]

Don't you see that you update Y as

Y(i,:) = y + h^2/6*(k0+2*k1+2*k2+k3);

instead of

Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);

when the ki were defined in the usual way as

k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);

?

[Jason Louison]

For solving differential equations that are not dependent on time explicity, you would do the following:

Say we have a differential equation dy/dt = f(y) = 1 + y. Since the variable dy/dt itself does not depend on time, but on y, you would do:

k1 = f(y(t))

k1 is just the differential equation evaluated at a given value of y. t is the time, but y is NOT dependent on time explicitly. It is to express the value of y at a given time t, y(t).

k2 = f(y(t)+0.5*h*k1)

k2 is the estimation of dy/dt = f(y) at a half of a time step forward. from f(y). Since k1 is the derivative of y, you can essentially approximate y(t+0.5h) by doing this. k3 takes it a step further by using the esitmated derivative k2 and repeating the process:

k3 = f(y(n)+0.5*h*k2)

k4 is the final iteration in this process. It estimates the value of dy/dt at one full time step away from f(y) by usin k3 as the derivate.

k4 = f(y+h*k3)

to get the approximation of y at the next time step, you use this formula:

y(n+1) = y(n)+ (h/6)*(k1*2k2*2k3+k4)

After obtaining this new value of y, the process is repeated. This is the method I use for numerically integrating functions that are not explicity a function of time. Here is an example code:

%Sprung Mass Damper Simulation%

k = 1000; %Spring Constant
b = 600; %Damping Constant
m = 1000; %Mass 
g = 9.81; %Gravity
p = k/m;
q = b/m;

tf = 120;
N = 10000;
h = tf/N;


t = linspace(0,tf,N+1);
x = zeros(1,N+1);
v = zeros(1,N+1);           
x(1) = 0;
v(1) = 0;
for n=1:N %fourth-order runge kutta method algorithm
    k1a = h*f(x(n),v(n),p,q,g);
    k1v = h*v(n);
    k2a = h*f(x(n)+1/2*k1v,v(n)+1/2*k1a,p,q,g);
    k2v = h*(v(n)+1/2*k1a);
    k3a = h*f(x(n)+1/2*k2v,v(n)+1/2*k2a,p,q,g);
    k3v = h*(v(n)+1/2*k2a);
    k4a = h*f(x(n)+k3v,v(n)+k3a,p,q,g);
    k4v = h*(v(n)+k3a);
    v(n+1) = v(n)+(1/6)*(k1a+2*(k2a+k3a)+k4a);
    x(n+1) = x(n)+(1/6)*(k1v+2*(k2v+k3v)+k4v);
end
a = -p*x-q*v-g;

% It was hard to see the displacement when all three functions were graphed
% on the same axis, since x is so mauch smaller than v & a
% I split them onto separate lines
subplot(3,1,1);
plot(t,x)
title('displacement')
subplot(3,1,2)
plot(t,v)
title('velocity')
subplot(3,1,3)
plot(t,a) %plot displacement, velocity, and acceleration vs time
title('acceleration')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = f(x,v,p,q,g) 

a = -p*x-q*v-g;
end

[Torsten]

Do you see the difference between the update of y in this code and the code this post is about ?

k1a = h*f(x(n),v(n),p,q,g);

k2a = h*f(x(n)+1/2*k1v,v(n)+1/2*k1a,p,q,g);

k3a = h*f(x(n)+1/2*k2v,v(n)+1/2*k2a,p,q,g);

k4a = h*f(x(n)+k3v,v(n)+k3a,p,q,g);

v(n+1) = v(n)+(1/6)*(k1a+2*(k2a+k3a)+k4a);

Here, you use 1/6 as multiplier.

k0 = h*f(t,y)';

k1 = h*f(t,y+k0*0.5)';

k2 = h*f(t,y+k1*0.5)';

k3 = h*f(t,y+k2)';

Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);

Here, you use h/6 as multiplier.

[Jason Louison]

Thank you for pointing out that mistype in the spring code, I was not aware of the error. Nevertheless these will give you the same result, it is just another way of writing the algorithm.

You can do it this way:

K1 = h*f(y(n)) K2 = h*f(y(n)+0.5*K1) K3 = h*f(y(n)+0.5*K2) K4 = h*f(y(n)+K3)

y(n+1) = (1/6)*(K1+2*(K2+K3)+K4)

OR, you can do it this way:

K1 = f(y(n)) K2 = f(y(n)+0.5*K1)) K3 = f(y(n)+0.5*K2)) K4 = f(y(n)+K3)

y(n+1) = (h/6)*(K1+2*(K2+K3)+K4)

I will correct the error in the code you have pointed out.

[Torsten]

K1 = f(y(n)) K2 = f(y(n)+0.5*h*K1)) K3 = f(y(n)+0.5*h*K2)) K4 = f(y(n)+h*K3)
y(n+1) = (h/6)*(K1+2*(K2+K3)+K4)

instead of

K1 = f(y(n)) K2 = f(y(n)+0.5*K1)) K3 = f(y(n)+0.5*K2)) K4 = f(y(n)+K3)
y(n+1) = (h/6)*(K1+2*(K2+K3)+K4)

It was a mistype in the chemical equilibrium code, not in the spring code.

And correcting it will lead to unphysical results.

[Jason Louison]

I have been busy for most of the day, but I will check out the error and correct it. Will update you via results soon.

[Jason Louison]

Torsten,

You are right, this is very strange. The strangest part is that I get normal, correct results when I write the algorithm as h^2/6*(k1+2*k2+2*k3+k4), but if I try writing the algorithm properly it, it breaks the simulation. It doesn't make any sense, as my other simulations work correctly with the algorithm I just explained (With the picture of the algorithm). Do you know what might be causing this? Attached below is the results of the "broken" simulation.

I will try re-setting the code back to the way you originally did it, but I am experienced in doing these simulations on spreadsheets, not code, and that is why I modified it the way I did. On a spreadsheet, using the runge kutta algorithm, you cannot include time in the alogorithm, even if a function changes with time. It would already factor it in. I did it this way in my spring damper system code as well, and it worked. Maybe because parameters are changing with time, I have to include time in the function. Anyways, I will get back to you soon with results of the new code.