# How can I make a 3D Plot with NDSolve?

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 == 10, Cx1 == 25,
Cx2 == 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 == 10, Cx1 == 25,
Cx2 == 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. >>

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 == 10, Cx1 == 25,
Cx2 == 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, Cx1, Cx2} /. solved[DL]], {DL, .01, .1}]

• I have tried to include the full code. My apologies if it is difficult to read. After attempting to take these steps I have received 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. >> – Nick L. Mar 15 '15 at 8:55
• I suggest you edit your code so it can easily be copied and pasted into Mathematica. You need to use SetDelayed and solved[DL] as shown in my edits. – djp Mar 15 '15 at 10:24
• Thank you very much. The code works very nicely. I also appreciate the tips which you gave me. I would up-vote you but apparently my reputation is too low. – Nick L. Mar 15 '15 at 18:52

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 == 10, Cx1 == 25,
Cx2 == 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] • Thank you for the code. I am very appreciative that you took the time to do this. This will be very helpful to me. – Nick L. Mar 15 '15 at 18:55