How can I plot these 3 ode45 functions in 1 graph?

Hello dear friends I have 3 codes which are different from each other with TH values and tspans ( for first TH=10 and tspan=[0 10], for second TH=20 & tspan=[0 20] and for the third one TH=40 & tspan=[0 40]). I need to plot all these 3 functions in one graph. At first step If it is possible, can you show me combining these three in one code and plot in same graph? Could you help me please? Codes are given below in order:

First:

X0=2;
Sin=1000;
So=0;
S0=So+Sin;
Ks=150;
Y=0.5;

XS0=[2, 1000];
tspan = [0 10]; 
[t,XS]=ode45(@BSfunction, tspan, XS0);
 
function dXSdt=BSfunction(TH,XS)
X = XS(1); S = XS(2);
X0=2;
Sin=1000;
So=0;
S0=So+Sin;
V=20;
TH=10;
ko=0.2;
kd=0.01;
Ks=150;
Y=0.5;
S=S0-(ko/Y)*X0*TH;
biomassgrowth=((ko*S*X)/(Ks+S))-(kd*X);
substratutilize = -((ko/Y)*X0);
dXSdt = [biomassgrowth; substratutilize];
end

Second:

X0=2;
Sin=1000;
So=0;
S0=So+Sin;
Ks=150;
Y=0.5;

XS0=[2, 1000];
tspan = [0 20]; 
[t,XS]=ode45(@BSfunction, tspan, XS0);
 
function dXSdt=BSfunction(TH,XS)
X = XS(1); S = XS(2);
X0=2;
Sin=1000;
So=0;
S0=So+Sin;
V=20;
TH=20;
ko=0.2;
kd=0.01;
Ks=150;
Y=0.5;
S=S0-(ko/Y)*X0*TH;
biomassgrowth=((ko*S*X)/(Ks+S))-(kd*X);
substratutilize = -((ko/Y)*X0);
dXSdt = [biomassgrowth; substratutilize];
end

Third:

X0=2;
Sin=1000;
So=0;
S0=So+Sin;
Ks=150;
Y=0.5;

XS0=[2, 1000];
tspan = [0 40]; 
[t,XS]=ode45(@BSfunction, tspan, XS0);
 
function dXSdt=BSfunction(TH,XS)
X = XS(1); S = XS(2);
X0=2;
Sin=1000;
So=0;
S0=So+Sin;
V=20;
TH=40;
ko=0.2;
kd=0.01;
Ks=150;
Y=0.5;
S=S0-(ko/Y)*X0*TH;
biomassgrowth=((ko*S*X)/(Ks+S))-(kd*X);
substratutilize = -((ko/Y)*X0);
dXSdt = [biomassgrowth; substratutilize];
end

Answer 1

[Star Strider]

The lines are not easily distinguishable.

Try this:

X0=2;
Sin=1000;
So=0;
S0=So+Sin;
Ks=150;
Y=0.5;
XS0=[2, 1000];

THv = [10 20 40];
for k = 1:numel(THv)
    tspan = [0 THv(k)];
    [t{k},XS{k}]=ode45(@(t,XS)BSfunction(t,XS,THv(k)), tspan, XS0);
end

figure
LineStl = {'-','-.','--'};
hold on
for k = 1:numel(THv)
    semilogy(t{k}, XS{k}, 'LineWidth',1.5, 'DisplayName', ['TH = ',num2str(THv(k))], 'LineStyle',LineStl{k})
end
hold off
grid
% lgdstr = compose('TH = %d', THv);
legend

figure
for k = 1:numel(THv)
    subplot(numel(THv), 1, k)
    plot(t{k}, XS{k}, 'LineWidth',1.5)
    title(compose('TH = %d', THv(k)));
    grid
end


function dXSdt=BSfunction(t,XS,TH)
X = XS(1); S = XS(2);
X0=2;
Sin=1000;
So=0;
S0=So+Sin;
V=20;
% TH=10;
ko=0.2;
kd=0.01;
Ks=150;
Y=0.5;
S=S0-(ko/Y)*X0*TH;
biomassgrowth=((ko*S*X)/(Ks+S))-(kd*X);
substratutilize = -((ko/Y)*X0);
dXSdt = [biomassgrowth; substratutilize];
end

The first figure plots them all on one axes, the second figure plots them each as separate subplots.

[Carey n'eville]

Thank you so much!!!!! I guess this is enough for me

[Star Strider]

As always, my pleasure!