function [ dy] = dydt( t,y,par )
%UNTITLED10 Summary of this function goes here
% Detailed explanation goes here
K1D=par(1);
C4=par(2);
KD=par(3);
tetaD=par(4); %procurar valor real
T=par(5);
KDBO=par(6);
Vs3=par(7);
fD5=par(8);
H=par(9);
K2D=par(10);
KNO3=par(11);
C2=par(12);
C1=par(13);
Ka=par(14);
K12= par(15);
KNIT= par(16);
GPI=par(17);
K1R= par(18);
SOD= par(19);
tetaa=par(20);
teta12=par(21);
tetaS=par(22);
teta1R=par(23);
Cs=par(24);
teta2D= par(25);
PNH3= par(26);
dy(1)= K1D*C4-KD*(tetaD^(T-20))*(y(2)/(KDBO+y(2)))*y(1)-(Vs3*(1-fD5)/H)*y(1) - (5/112)*K2D*((teta2D)^(T-20))*(KNO3/(KNO3+y(2)))*C2;
dy(2)= Ka*(tetaa^(T-20))*(Cs-y(2))-(KD*(tetaD^(T-20)))*(y(2)/(KDBO+ y(2)))*y(1)- (64/14)*K12*(teta12^(T-20))*(y(2)/(KNIT+ y(2)))*C1-(32/12)*K1R*(teta1R^(T-20))*C4+(GPI*((32/12)+(48/12)*(1-PNH3)))*C4 - (SOD/H)*(tetaS^(T-20));
dy=[dy(1) dy(2)];
end
close all
clear all
clc
K1D=0.02;
C4=0.0;
KD=0.15;
tetaD=1.047;
T=25;
KDBO=0.50;
Vs3=25920.0;
fD5=0.5;
H=725.0;
K2D=0.09;
KNO3=0.10;
C2=2.0;
C1=0.50;
Ka=0.7;
K12=0.02;
KNIT=0.3;
GPI=2.30;
K1R=0.16;
SOD=0.5;
tetaa=1.024;
teta12=1.012;
tetaS=1.08;
teta1R=0.0;%não tem
Cs=0.0; %no tempo inicial é zero
teta2D=1.045;
PNH3=5000;
t = 0.0 : 0.150 : 10.0;
y0=[0 8.33];
par=[K1D C4 KD tetaD T KDBO Vs3 fD5 H K2D KNO3 C2 C1 Ka K12 KNIT GPI K1R SOD tetaa teta12 tetaS teta1R Cs teta2D PNH3];
modelo= @(t,y) dydt(t,y,par);
[t, y]= ode45(modelo,t,y0);
figure(1)
plot(t,y(:,1), t,y(:,2));
The program rotates, but I get negative values, which can not occur since the variable under study is the concentration of oxygen.