How can I make a 3D Plot with NDSolve?

Question

I have a system of coupled ODEs which I am solving using Mathematica. I am solving these using NDSolve:

solved=NDSolve[{Cs'[t] == DL*Stuff, 
  Cx1'[t] == DL*Stuff, 
  Cx2'[t] == DL*Stuff, Cs[0] == 10, Cx1[0] == 25, 
  Cx2[0] == 7}, {Cs, Cx1, Cx2}, {t, 0, 1000}]

Where DL is a constant. Generally following this I will plot them as follows:

Plot[Evaluate[{Cs[t], Cx1[t], Cx2[t]} /. solved], {t, 0, 1000}, 
 PlotRange -> All]

However, what I would like to do is make a 3D plot which shows how the solution to the system changes for different values of the constant DL. How can I do this?

Constraint: The equations mentioned above cannot be solved analytically, or decoupled or whatever.

Full code after attempting changes:

(First, I will Define Constants)

Subscript[Y, sx1] = 1/.14;
U1max = .5;
Subscript[KM, 1] = Subscript[KM, 2] = 10; 
Subscript[Y, x1x2] = 2 ;
U2max = .11;

(Next, I will define the prior algebraic expressions)

Clear[Cs]

U1 = (U1max*Cs[t])/(Subscript[KM, 1] + Cs[t]);

U2 = (U2max*Cs[t])/(Subscript[KM, 2] + Cs[t]);

(Now I put in the remaining algebraic expressions)

rgx1 = U1*Cx1[t];

rgx2 = U2*Cx2[t];

(Now we can enter the coupled ODEs)

solved[DL_] = 
 NDSolve[{Cs'[t] == DL*(250 - Cs[t]) - Subscript[Y, sx1]*rgx1, 
   Cx1'[t] == -DL*Cx1[t] + rgx1 - Subscript[Y, x1x2]*rgx2, 
   Cx2'[t] == -DL*Cx2[t] + rgx2, Cs[0] == 10, Cx1[0] == 25, 
   Cx2[0] == 7}, {Cs, Cx1, Cx2}, {t, 0, 1000}]

Plot3D[Evaluate[Cx2[t] /. solved], {t, 0, 1}, {DL, .01, .1}, 
 PlotRange -> All]

I get two errors:

(1) NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

(2) ReplaceAll::reps: {solved} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

Answer 1

Here you go:

solved[DL_] := 
 NDSolve[{Cs'[t] == DL*(250 - Cs[t]) - Subscript[Y, sx1]*rgx1, 
  Cx1'[t] == -DL*Cx1[t] + rgx1 - Subscript[Y, x1x2]*rgx2, 
  Cx2'[t] == -DL*Cx2[t] + rgx2, Cs[0] == 10, Cx1[0] == 25, 
  Cx2[0] == 7}, {Cs, Cx1, Cx2}, {t, 0, 1000}]

Plot3D[Cx2[tt] /. solved[DL]], {tt, 0, 1000}, {DL, 0.01, 0.1}, 
  PlotRange -> All]

Some important points or it won't work:

(1) the value tt for time has to have a different name to the value t for time inside the NDSolve

(2) the SetDelayed in the definition of solved is important

(3) it's horribly inefficient

But it works nicely.

Comment:

Usually the above approach isn't very useful; e.g. you would find it more useful to plot the final state as a function of DL.

Plot[Evaluate[{Cs[1000], Cx1[1000], Cx2[1000]} /. solved[DL]], {DL, .01, .1}]

Answer 2

Rewriting the code:

Ysx1 = 1/.14;
U1max = .5;
KM1 = KM2 = 10;
Yx1x2 = 2;
U2max = .11;
Clear[Cs]
U1 = (U1max*Cs[t])/(KM1 + Cs[t]);
U2 = (U2max*Cs[t])/(KM2 + Cs[t]);
rgx1 = U1*Cx1[t];
rgx2 = U2*Cx2[t];
sol = ParametricNDSolve[{Cs'[t] == DL*(250 - Cs[t]) - Ysx1*rgx1, 
   Cx1'[t] == -DL*Cx1[t] + rgx1 - Yx1x2*rgx2, 
   Cx2'[t] == -DL*Cx2[t] + rgx2, Cs[0] == 10, Cx1[0] == 25, 
   Cx2[0] == 7}, {Cs, Cx1, Cx2}, {t, 0, 1000}, {DL}]

Assuming the desired 3D plot relates to {DL, t, Cx2}:>

f[DL_, t_] := (Cx2[DL] /. sol)[t]
Plot3D[Evaluate[f[x, y]], {y, 0, 1000}, {x, 0.01, 0.1}, 
PlotRange -> All, AxesLabel -> {"t", "DL", "Cx2"}, Mesh -> False]