Trouble solving ODE equations

I am trying to solve these differential equations in order to plot graphs, however, when running the code I am receiving multiple errors but I'm unsure why as I can't see why the code is not running smoothly? I'm looking for some insight on why this might be happening and how I can solve the problem?

Equations and Parameter values:

function dydt = Dunster(t,y)

%Parameters and their values%

K1a = 20;

Y1a = 1;

K1b = 0.25;

K2a = 14;

K2am = 72;

K2b = 2.6;

K2bm = 10;

K2c = 24;

K2cm = 20;

K3a = 10;

K3b = 0.05;

K3c = 24;

K3cm = 20;

K4a = 2.3;

K4am = 58;

K4b = 2000;

K4bm = 210;

K4c = 1.3;

K5a = 0.0014;

K5b = 0.35;

K6 = 2000;

K6m = 2000;

Kx = 50;

Kb = 0.01;

Lb = 1;

Lx = 1;

Bxva = 0;

H = K1a*Y1a*exp(-Y1a*t);

bl = 0.5*((Kb+Lb+Bxva)-sqrt((Kb+Lb+Bxva)^2 - 4*Lb*Bxva));

%Equation 1%

%Generation and inactivation of factor Xa%

F1= y(1);

%y(1)= -K6*F1;

% Factor I

%Equation 2%

F5=y(2);

%y(2)=(-K2a*F2a*F5)/(F5+K2am(1+F1/K6m) - ((K2b*F10a*F5)/F5+K2bm(1+F2/K4am)));

%Equation 3%

F2=y(3);

%y(3)= ((-K4a*F1*F2)/F2+K4am(1+F5/K2bm))-((K4b*b1*F2)/(F2+K4bm));

%Equation 4%

C=y(4);

%y(4)= (-K5a*C);

%Equation 5%

F10a=y(5);

%y(5)= H +(K3c*Ca*Bxva)/(Bxva+K3cm) - K1b*F10a - K3a*F10a*F5a;

%Equation 6%

F5a=y(6);

%y(6)= (K2a*F2a*F5)/(F5+K2am(1+F1/K6m)) + (K2b*F10a*F5a)/(F5+K2bm(1+F1/K4am))+ K3b*Bxva - (K2c*Ca*F5a)/(F5a+K2am) - K3a*F10a*F5a;

%Equation 7%

Bxva=y(7);

%y(7)= K3a*F10a*F5a - K3b*Bxva - (K3c*Ca*Bxva)/(Bxva+K3cm);

%Equation 8%

F2a=y(8);

%y(8)= (K4a*F1*F2)/(F2+K4am(1+F5/K2bm)) + (K4b*BF2)/(F2+K4bm) - K4c*F2a;

%Equation 9%

Ca=y(9);

%y(9)= K5a*C - K5b*Ca;

%Equation 10%

F10i=y(10);

%y(10)= K1b*F10a + K3b*Bxva;

%Equation 11%

F5i=y(11);

%y(11) = (K2c*Ca*F5a)/(F5a+K2cm) + (K3c*Ca*Bxva/(Bxva+K3cm));

%Equation 12%

F2i=y(12);

%y(12) = K4c*F2a;

%Equation 13%

Ci=y(13);

%y(13) = K5b*Ca;

%Equation 14%

F1a=y(14);

%y(14) = K6*F1;

dydt=[

-K6*F1; % d F1 / dt

((-K2a*F2a*F5)/(F5+K2am*(1+F1/K6m)) - ((K2b*F10a*F5)/F5+K2bm*(1+F2/K4am))); % d F5 / dt

((-K4a*F1*F2)/(F2+K4am*(1+F5/K2bm))-((K4b*bl*F2)/(F2+K4bm))); % d F2 / dt

-K5a*C; % d C / dt

H+(K3c*Ca*Bxva)/(Bxva+K3cm) - K1b*F10a - K3a*F10a*F5a ;

% d F10a / dt

(K2a*F2a*F5)/(F5+K2am(1+F1/K6m)) + (K2b*F10a*F5a)/(F5+K2bm(1+F1/K4am))

+ K3b*Bxva - (K2c*Ca*F5a)/(F5a+K2am) - K3a*F10a*F5a; % d F5a / dt

K3a*F10a*F5a - K3b*Bxva - (K3c*Ca*Bxva)/(Bxva+K3cm); % d Bxva / dt

(K4a*F1*F2)/(F2+K4am(1+F5/K2bm)) + (K4b*BF2)/(F2+K4bm) - K4c*F2a ;

% d F2a / dt

K5a*C - K5b*Ca ; % d Ca / dt

K1b*F10a + K3b*Bxva ; % d F10i / dt

((K2c*Ca*F5a)/(F5a+K2cm) + (K3c*Ca*Bxva/(Bxva+K3cm))) ; % d F5i / dt

K4c*F2a ; % d F2i / dt

K5b*Ca ; % d Ci / dt

K6*F1 ; % d F1a / dt

];

end

Running Part of the Code:

tspan=[0 60];

F1=10000;

F5=30;

F2=1000;

C=100;

F10a=0;

F5a=0;

Bxva=0;

F2a=0;

Ca=0;

F10i=0;

F5i=0;

F2i=0;

Ci=0;

F1a=0;

y0=[F1;F5;F2;C;F10a;F5a;Bxva;F2a;Ca;F10i;F5i;F2i;Ci;F1a];

[t,y]=ode45(@Dunster, tspan, y0);%,[], pars);

plot(t/60, y(:,14)*1000, "r")

xlabel("Time [Min]")

ylabel("Concentration [nM]")

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[Star Strider]

There were several missing multiplication operators, and a reference to ‘BF2’ that I corrected to ‘F2’, since there is no ‘B’ that I can find, so no missing multiplication operator there. With those edits, and concatenating an additional 0 to ‘y0’ to make its length equal to the number of differential equations, this runs without error:

function dydt = Dunster(t,y)
%Parameters and their values%
K1a = 20;
Y1a = 1;
K1b = 0.25;
K2a = 14;
K2am = 72;
K2b = 2.6;
K2bm = 10;
K2c = 24;
K2cm = 20;
K3a = 10;
K3b = 0.05;
K3c = 24;
K3cm = 20;
K4a = 2.3;
K4am = 58;
K4b = 2000;
K4bm = 210;
K4c = 1.3;
K5a = 0.0014;
K5b = 0.35;
K6 = 2000;
K6m = 2000;
Kx = 50;
Kb = 0.01;
Lb = 1;
Lx = 1;
Bxva = 0;
H = K1a*Y1a*exp(-Y1a*t);
bl = 0.5*((Kb+Lb+Bxva)-sqrt((Kb+Lb+Bxva)^2 - 4*Lb*Bxva));
%Equation 1%
%Generation and inactivation of factor Xa%
F1= y(1);
%y(1)= -K6*F1;
% Factor I
%Equation 2%
F5=y(2);
%y(2)=(-K2a*F2a*F5)/(F5+K2am(1+F1/K6m) - ((K2b*F10a*F5)/F5+K2bm(1+F2/K4am)));
%Equation 3%
F2=y(3);
%y(3)= ((-K4a*F1*F2)/F2+K4am(1+F5/K2bm))-((K4b*b1*F2)/(F2+K4bm));
%Equation 4%
C=y(4);
%y(4)= (-K5a*C);
%Equation 5%
F10a=y(5);
%y(5)= H +(K3c*Ca*Bxva)/(Bxva+K3cm) - K1b*F10a - K3a*F10a*F5a;
%Equation 6%
F5a=y(6);
%y(6)= (K2a*F2a*F5)/(F5+K2am(1+F1/K6m)) + (K2b*F10a*F5a)/(F5+K2bm(1+F1/K4am))+ K3b*Bxva - (K2c*Ca*F5a)/(F5a+K2am) - K3a*F10a*F5a;
%Equation 7%
Bxva=y(7);
%y(7)= K3a*F10a*F5a - K3b*Bxva - (K3c*Ca*Bxva)/(Bxva+K3cm);
%Equation 8%
F2a=y(8);
%y(8)= (K4a*F1*F2)/(F2+K4am(1+F5/K2bm)) + (K4b*BF2)/(F2+K4bm) - K4c*F2a;
%Equation 9%
Ca=y(9);
%y(9)= K5a*C - K5b*Ca;
%Equation 10%
F10i=y(10);
%y(10)= K1b*F10a + K3b*Bxva;
%Equation 11%
F5i=y(11);
%y(11) = (K2c*Ca*F5a)/(F5a+K2cm) + (K3c*Ca*Bxva/(Bxva+K3cm));
%Equation 12%
F2i=y(12);
%y(12) = K4c*F2a;
%Equation 13%
Ci=y(13);
%y(13) = K5b*Ca;
%Equation 14%
F1a=y(14);
%y(14) = K6*F1;
dydt=[-K6*F1; % d F1 / dt
    
    ((-K2a*F2a*F5)/(F5+K2am*(1+F1/K6m)) - ((K2b*F10a*F5)/F5+K2bm*(1+F2/K4am))); % d F5 / dt
    
    ((-K4a*F1*F2)/(F2+K4am*(1+F5/K2bm))-((K4b*bl*F2)/(F2+K4bm))); % d F2 / dt
    
    -K5a*C; % d C / dt
    
    H+(K3c*Ca*Bxva)/(Bxva+K3cm) - K1b*F10a - K3a*F10a*F5a ; 
    % d F10a / dt
    
    (K2a*F2a*F5)/(F5+K2am*(1+F1/K6m)) + (K2b*F10a*F5a)/(F5+K2bm*(1+F1/K4am))
+ K3b*Bxva - (K2c*Ca*F5a)/(F5a+K2am) - K3a*F10a*F5a; % d F5a / dt 
K3a*F10a*F5a - K3b*Bxva - (K3c*Ca*Bxva)/(Bxva+K3cm); % d Bxva / dt
(K4a*F1*F2)/(F2+K4am*(1+F5/K2bm)) + (K4b*F2)/(F2+K4bm) - K4c*F2a ;
% d F2a / dt
K5a*C - K5b*Ca ; % d Ca / dt
K1b*F10a + K3b*Bxva ; % d F10i / dt
((K2c*Ca*F5a)/(F5a+K2cm) + (K3c*Ca*Bxva/(Bxva+K3cm))) ; % d F5i / dt
K4c*F2a ; % d F2i / dt
K5b*Ca ; % d Ci / dt
K6*F1 ; % d F1a / dt
];
end

% % Running Part of the Code:
tspan=[0 60];
F1=10000;
F5=30;
F2=1000;
C=100;
F10a=0;
F5a=0;
Bxva=0;
F2a=0;
Ca=0;
F10i=0;
F5i=0;
F2i=0;
Ci=0;
F1a=0;
y0=[F1;F5;F2;C;F10a;F5a;Bxva;F2a;Ca;F10i;F5i;F2i;Ci;F1a;0];
[t,y]=ode45(@Dunster, tspan, y0);%,[], pars);
plot(t/60, y(:,14)*1000, "r")
xlabel("Time [Min]")
ylabel("Concentration [nM]")

.

[Pr0t0nZ]

Thank you for your reply, I have now got the code running but is only producing 1 graph instead of 14 (1 for each equation), is there any chance you know how to get MATLAB to produce a graph for each equation?

Really appreciate your help with this, Thank you.

[Star Strider]

My pleasure!

Put this after the ode45 call:

figure
for k = 1:size(y,2)
    subplot(8,2,k)
    plot(t/60, y(:,k)*1000, "r")
    xlabel("Time [Min]")
    ylabel("Concentration [nM]")
    title(sprintf('y_{%d}',k))
end
pos = get(gcf, 'Position');
set(gcf, 'Position',pos+[0 -500 0 500])

The rest of the code is unchanged.

There are other options (such as plotting them all on the same axes), however with 15 variables with significantly different magnitudes, that would be difficult to interpret.

Experiment to get the result you want.

[Pr0t0nZ]

Hi,
Thank you for this insight, i got multiple graphs (15) but all of the graphs are sharing the window so they are crushed up and small. How can I get the full size of each graph?

Cheers,

Jack

[Star Strider]

My pleasure!

Fro the plots:

for k = 1:size(y,2)
    figure
    plot(t/60, y(:,k)*1000, "r")
    xlabel("Time [Min]")
    ylabel("Concentration [nM]")
    title(sprintf('y_{%d}',k))
    set(gca, 'XScale','log')                                        % Optional
end

Also consider:

y0=[F1;F5;F2;C;F10a;F5a;Bxva;F2a;Ca;F10i;F5i;F2i;Ci;F1a;0]+eps;     % Adding ‘eps’ Optional, Prevents Some Variables From Being Identically 0

.

[Pr0t0nZ]

Hi,
I just wanted to say a massive thank you for all your help you have given me, honestly I can't stress enough how long i've been trying to sort this problem with no results and now my code is fully working and producing graphs!
So once again, thank you so much for all your help!
Cheers,
Jack

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.