# Solving and plotting ODEs while varying one of the initial conditions

How can I vary an initial condition in the numerical solution of a system of ODEs and then make a 3D plot of the solution space with that condition as one of the variables.

sol = NDSolve[
{y''[t] + y'[t] + 4*y[t] + x[t] - x'[t] == 0,
x'[t] + 3*c*y[t] == 0,
y == 1, y' == 1,
x == 1}, {x, y}, {t, 0, 10}]


I want to vary one of the condition variables, say y over u = Range[-10, 10, step]. Then I want to make a 3D plot of the solution space (x[t[, y[t], y[u]).

Can anyone please guide me in solving this query.

• Your equation can be solved analytically, so you can just exchange y==1 with y==z and DSolve the equation and Plot3D[x[t] /. sol, {t, -1, 1}, {z, -1, 1}]. Mar 7, 2013 at 8:09
• @xzczd In fact, the equation in my problem involves multipliers in terms of time, I have put the simplest form. I need a way to solve numerically. Mar 7, 2013 at 8:44
• @MuhammadZubair, have a look at ParametricNDSolve and see if that helps.
– user21
Mar 7, 2013 at 9:00
• @ruebenko Wow, I think I'd better upgrade to v9 quickly. Mar 7, 2013 at 9:06
• @xzczd, if you need a parametric version of NDSolve this is it; automatic sensitivity computation. It's pretty cool, I think. Also a much better Event language and better DAE solving capabilities, are just a few highlights of V9 NDSolve. Hope that's enough of a teaser...
– user21
Mar 7, 2013 at 10:20

Playing with @ruebenko's suggestion,

c = 1
pf = ParametricNDSolveValue[
{y''[t] + y'[t] + 4*y[t] + x[t] - x'[t] == 0,
x'[t] + 3*c*y[t] == 0, y == 1,
y' == u,
x == 1}, y, {t, 0, 10}, {u \[Element] Reals}]

Plot3D[pf[u][t], {u, -10, 10}, {t, 0, 10}, PlotRange -> All] ruebenko has pointed out the best way to solve this in the comment above: use ParametricNDSolve in version 9. Then, I'd like to post my clumsy solution with version 8 since I've already finished it…:

c = 1;
zmin = -1; zmax = 1; n = 25;
tmin = -1; tmax = 1;
ex =
Table[Join[{z}, #] & /@
Transpose[{First[(x /. First@#)["Coordinates"]], (x /. First@#)["ValuesOnGrid"]}] &@
NDSolve[{y''[t] + y'[t] + 4 y[t] + x[t] - x'[t] == 0,
x'[t] + 3 c y[t] == 0, y == z, y' == 1, x == 1},
{x, y}, {t, tmin, tmax}],
{z, zmin, zmax, (zmax - zmin)/n}];

ListPlot3D[Flatten[ex, 1]] Just for your sample, it can be solved analytically, too:

c = 1;
sol = DSolve[{y''[t] + y'[t] + 4 y[t] + x[t] - x'[t] == 0,
x'[t] + 3 c y[t] == 0, y == z, y' == 1, x == 1},
{x, y}, {t}];
Plot3D[x[t] /. sol, {z, -1, 1}, {t, -1, 1}] 