How can i run this script to make a plot for a lot of T

function [f] = acidpretreatment(t,x,T,xyl_zero)

%Components

xyl= x(1);

arab= x(2);

acet= x(3);

acet_dot = x(4);

cel= x(5);

Xylo= x(6);

Acetic= x(7);

Acetic_dot = x(8);

Arabi = x(9);

glu = x(10);

furf= x(11);

%---------------------------------------------------------------------

%Data

R=8.314;

HA=0.07; % m=948.320,38 kg

i=100:180

T=375+1i; % kelvin "Pressure 10 bar"

k10=2.37; k20=2.17; k1ac=2.37; k30=2.37; %min^-1%

n1=1.51; n2=0.29;n1ac=0.604;n3=1.359;

E1=83.3*1e3;E2=143.5*1e3;E1ac=83.3*1e3;E3=94.962*1e3; %J/mol

Bmax=120;

Xc=0.40680; Xh1=0.22136; Xh2=0.03786;Xh3=0.03330; %wt

xylmax=Xh1*Bmax; arabmax=Xh2*Bmax;acetmax=Xh3*Bmax;celmax=Xc*Bmax; %g/L

%---------------------------------------------------------------------

% Effective coefficients

ef1=1/(1+(xyl/xylmax)^20);

ef2=1/(1+(arab/arabmax)^20);

ef3=1/(1+(acet/acetmax)^20);

ef4=1/(1+(cel/celmax)^20);

% Effective concentrations

xylef=0.95*ef1*xyl+((1-ef1)*xylmax);

arabef=0.95*ef2*arab+((1-ef2)*arabmax);

acetef=0.95*ef3*acet+((1-ef3)*acetmax);

celef=0.95*ef4*cel+((1-ef4)*celmax);

% K min^(-1)

K1=(k10*1e10)*(HA^n1)*exp(-E1/(R*T));

K2=(k20*1e15)*(HA^n2)*exp(-E2/(R*T));

K3=(k30*1e10)*(HA^n3)*exp(-E3/(R*T)) ;

K1ac=k1ac*(10^10)*(HA^n1ac)*(xyl/xyl_zero)^2*exp(-E1ac/(R*T));

% Rates

rxyl=-K1*xylef;

rarab=-K1*arabef;

racet=-K1ac*acetef;

rcel=-K3*celef;

rXylo=(1/0.88)*K1*xylef-K2*Xylo;

rArabi=(1/0.88)*K1*arabef-K2*Arabi;

rAcetic=(1/0.7)*K1ac*acetef;

rglu=(1/0.9)*K3*celef;

% Mass balance

dxyldt=rxyl;

darabdt=rarab;

dacetdt = acet_dot;

d2acetdt2 = (racet - 10*acet_dot)/400;

daceticdt = Acetic_dot;

d2aceticdt2 = (rAcetic - 10*Acetic_dot)/400;

dceldt=rcel;

dXylodt=rXylo;

dArabidt=rArabi;

dgludt=rglu;

dfurfdt=0.64*K2*(Xylo+Arabi);

%-----------------------------------

f =[ dxyldt; darabdt; dacetdt; d2acetdt2; dceldt; dXylodt; daceticdt; d2aceticdt2; dArabidt; dgludt; dfurfdt];

end

The Xo is another file that i run first and then I run then on cmd : [t, x] = ode45(@(t,y) acidpretreatment(t,y,xyl_zero), [0 240], [x0]);

How it is possible to sovle this ode for T from 100 to 180.

Answer 1

[VBBV]

T = 100:10:180;
x0 =100.8;
xyl_zero = 4;
% x0 = rand(1,length(T))
for i = 1:length(T)
[t,x] = ode45(@acidpretreatment, [0 24], x0);
plot(t,x)
hold on
end 

function [dxyldt darabdt dacetdt d2acetdt2 daceticdt d2aceticdt2 dceldt dXylodt dArabidt dgludt dfurfdt] = acidpretreatment(T,xyl_zero)
%Components
X = rand(1,11);
xyl=    X(1);
arab=   X(2);
acet=   X(3);
acet_dot = X(4);
cel=    X(5);
Xylo=   X(6);
Acetic= X(7);
Acetic_dot = X(8);
Arabi = X(9);
glu = X(10);
furf=   X(11);
 %---------------------------------------------------------------------
 %Data 
 R=8.314;     
 HA=0.07;   % m=948.320,38 kg
%i=100:180  
 T=375+T; % kelvin "Pressure 10 bar"  
 k10=2.37;
 k20=2.17; 
 k1ac=2.37; 
 k30=2.37; %min^-1%
 n1=1.51;
 n2=0.29;
 n1ac=0.604;
 n3=1.359;
 E1=83.3*1e3;E2=143.5*1e3;E1ac=83.3*1e3;E3=94.962*1e3; %J/mol   
 Bmax=120; 
 
Xc=0.40680; Xh1=0.22136; Xh2=0.03786;Xh3=0.03330; %wt
xylmax=Xh1*Bmax; arabmax=Xh2*Bmax;acetmax=Xh3*Bmax;celmax=Xc*Bmax; %g/L
%---------------------------------------------------------------------
%  Effective coefficients
ef1=1/(1+(xyl/xylmax)^20);
ef2=1/(1+(arab/arabmax)^20);
ef3=1/(1+(acet/acetmax)^20);
ef4=1/(1+(cel/celmax)^20);
%  Effective concentrations
xylef=0.95*ef1*xyl+((1-ef1)*xylmax);
arabef=0.95*ef2*arab+((1-ef2)*arabmax);
acetef=0.95*ef3*acet+((1-ef3)*acetmax);
celef=0.95*ef4*cel+((1-ef4)*celmax);
% K min^(-1)
K1=(k10*1e10)*(HA^n1)*exp(-E1/(R*T)); 
K2=(k20*1e15)*(HA^n2)*exp(-E2/(R*T)); 
K3=(k30*1e10)*(HA^n3)*exp(-E3/(R*T)) ; 
K1ac=k1ac*(10^10)*(HA^n1ac)*(xyl/xyl_zero)^2*exp(-E1ac./(R*T));
% Rates

rxyl=-K1*xylef;
rarab=-K1*arabef;
racet=-K1ac*acetef;
rcel=-K3*celef;
rXylo=(1/0.88)*K1*xylef-K2*Xylo;
rArabi=(1/0.88)*K1*arabef-K2*Arabi;
rAcetic=(1/0.7)*K1ac*acetef;
rglu=(1/0.9)*K3*celef;
% Mass balance

dxyldt=rxyl;
darabdt=rarab;
dacetdt = acet_dot;
d2acetdt2 = (racet - 10*acet_dot)/400;
daceticdt = Acetic_dot;
d2aceticdt2 = (rAcetic - 10*Acetic_dot)/400;
dceldt=rcel;
dXylodt=rXylo;
dArabidt=rArabi;
dgludt=rglu;
dfurfdt=0.64*K2*(Xylo+Arabi);
 
%-----------------------------------
%f = [dxyldt; darabdt; dacetdt; d2acetdt2; dceldt; dXylodt; daceticdt; d2aceticdt2; dArabidt; dgludt; dfurfdt]
end

Change the number of return arguments to function and try it.