MATLAB: Can someone check the code and rectify me, getting constant output but needs changing.

differential equationserrorfunctionsMATLABodesolve

Can someone check my code and rectify my error.
Here I have given a main code, two functions. On running it gives me constant output but I need that value of 'a' in main function should keep changing with each loop and final value I need like Glc get stored, at end of loop, to plot a variable graph. I have tried to simplify code and added possible change in comments format, but there's still problem. So, I needed help in it. In advance Thanks a lot for helping me out.
clear
close all
clc
M = {'Gal','Urd','Glc','Asn','Amm','Glu','Lac','Gln','Asp'};
v = 100;
cellden = 200000;
mb = 2.76709;
umax = 0.0650;
gal = 0;
urd = 0;
glc = 36;
asn = 6.09;
amm = 0.1;
glu = 2.125333333;
lac = 1.391166667;
gln = 0.01;
asp = 1.05;
for k = 1:216
a = flim(glc,asn);
b = 2;
u = umax*a*b;
udeathf = 3;
qAmm = 12;
qUrd = 19;
qGln = 15;
qGlc = -10;
qGal = 11;
qLac = 12;
qAsn = -13;
qAsp = -12;
qGlu = -9;;
tspan = [0 216];
xo = [100;20000000;0;0;3600;609;10;212.533;139.11667;1;105;276.709];
[t,X]=ode23s(@(t,X) new_try_func(t,X,v,cellden,mb,u,udeathf,gal,urd,glc,asn,amm,glu,lac,gln,asp,qAmm,qUrd,qGln,qGlc,qGal,qLac,qAsn,qAsp,qGlu),tspan,xo);
cd(k) = X(1,2)/X(1,1);
time(k) = k;
v = X(1,1);
cellden = X(1,2)/X(1,1);
gal = X(1,3)/X(1,1);
urd = X(1,4)/X(1,1);
glc = X(1,5)/X(1,1);
asn = X(1,6)/X(1,1);
amm = X(1,7)/X(1,1);
glu = X(1,8)/X(1,1);
lac = X(1,9)/X(1,1);
gln = X(1,10)/X(1,1);
asp = X(1,11)/X(1,1);
mb = X(1,12)/X(1,1);
end
%cell_density(:,1) = X(:,2)./X(:,1);
%Gal(:,1) = X(:,3)./X(:,1);
%Urd(:,1) = X(:,4)./X(:,1);
%Glc(:,1) = X(:,5)./X(:,1);
%Asn(:,1) = X(:,6)./X(:,1);
%Amm(:,1) = X(:,7)./X(:,1);
%Glu(:,1) = X(:,8)./X(:,1);
%Lac(:,1) = X(:,9)./X(:,1);
%Gln(:,1) = X(:,10)./X(:,1);
%Asp(:,1) = X(:,11)./X(:,1);
%mAb(:,1) = X(:,12)./X(:,1);
plot(time,cd); % or plot(cell_density);
xlim([0 216]);
Next are function code in simplified way(not original)
function dydt = new_try_func(~,X,v,cellden,mb,u,udeathf,gal,urd,glc,asn,amm,glu,lac,gln,asp,qAmm,qUrd,qGln,qGlc,qGal,qLac,qAsn,qAsp,qGlu)
mmAb = 0.410;
YmAb_X = 3.39;
dydt = zeros(12,1);
Fin = 0.1;
Fout = 0;
Gal = 0;
Urd = 0;
Glc = 144.37;
Asn = 26.99;
Amm = 0.06;
Glu = 12.19;
Lac = 16.64;
Gln = 0;
Asp = 51.95;
dydt(1,1) = Fin - Fout; % volume
dydt(2,1) = ((u - udeathf)*v*cellden) - (Fout*cellden); %cell density
dydt(3,1) = (Fin*Gal) - (Fout*gal) + (qGal*v*cellden); % Gal
dydt(4,1) = (Fin*Urd) - (Fout*urd) + (qUrd*v*cellden); % Urd
dydt(5,1) = (Fin*Glc) - (Fout*glc) + (qGlc*v*cellden); % Glc
dydt(6,1) = (Fin*Asn) - (Fout*asn) + (qAsn*v*cellden); % Asn
dydt(7,1) = (Fin*Amm) - (Fout*amm) + (qAmm*v*cellden); % Amm
dydt(8,1) = (Fin*Glu) - (Fout*glu) + (qGlu*v*cellden); % Glu
dydt(9,1) = (Fin*Lac) - (Fout*lac) + (qLac*v*cellden); % Lac
dydt(10,1) = (Fin*Gln) - (Fout*gln) + (qGln*v*cellden); % Gln
dydt(11,1) = (Fin*Asp) - (Fout*asp) + (qAsp*v*cellden); % Asp
dydt(12,1) = (-1*Fout*mb) + ((YmAb_X*u + mmAb)*v*cellden); %mAb
%v = X(1);
%cellden = X(2)/X(1);
%gal = X(3)/X(1);
%urd = X(4)/X(1);
%glc = X(5)/X(1);
%asn = X(6)/X(1);
%amm = X(7)/X(1);
%glu = X(8)/X(1);
%lac = X(9)/X(1);
%gln = X(10)/X(1);
%asp = X(11)/X(1);
%mb = X(12)/X(1);
end
Function for 'a'
function a = flim(glc,asn)
KGlc = 14;
KAsn = 2;
a = ((glc*asn)/((glc + KGlc)*(asn + KAsn)));
end

Best Answer

[t,X]=ode23s(@(t,X) new_try_func(t,X,v,cellden,mb,u,udeathf,gal,urd,glc,asn,amm,glu,lac,gln,asp,qAmm,qUrd,qGln,qGlc,qGal,qLac,qAsn,qAsp,qGlu),tspan,xo);
There is one row of X for each t value. The rows correspond to the X vector at different time steps.
glc = X(1,5)/X(1,1);
The different rows correspond to the X vector at different time steps. The first row, X(1,:), corresponds to the X values at the first time step. Which is to say, it is always going to be the same as the boundary conditions you pass in.
You are therefore going to be setting the new value of glc the same as what you passed in for the boundary conditions. You are passing in 100 for x0(1) and 3600 for x0(6) so you will get out 3600/100 = 36 which is the same as your initial glc. Nothing can change.
You should be considering extracting from (for example) X(end,5)/X(end,1)
(My internal tracking number for this is T0098912)
Related Question