# NDSolve with two parameters

I was trying to solve a ODE numerically. It has two parameters (w and z0) which I want to vary. The following code gives an error: "Cannot find starting value for the variable".

g[z_] := 1 + 3 ( (z/z0))^4 - 4 ( (z/z0))^3
DE = x''[z] + D[g[z]/z^2, z]/(g[z]/z^2) x'[z] +  (w^2)/(g[z])^2 x[z] == 0;

a = 10^-4;
zb = 0.0001;

pfun = ParametricNDSolveValue[{
DE, x'[z0 - a] == - I (z0 w)/(6 w^2) x[z0 - a], x[zb] == 1
}, x'[zb], {z, zb, z0 - a}, {w, z0}]
Plot[Evaluate[Table[
Re[pfun[w, z0]]
, {z0, 0.5, 1, .1}]]
, {w,-0.01,0.01}, PlotRange -> All]

• This error is an "Initial Condition" definition error. Please consider rewriting x'[z0 - a] == - I (z0 w)/(6 w^2) x[z0 - a], x[zb] == 1 in a more clear way. Nov 10, 2015 at 7:16
• @SerhanAya Thanks for your comment! But I have no idea how to modify it. Can you suggest something which might work? Nov 10, 2015 at 11:09

This seems to be caused by the boundary condition evaluated at a variable point. As an illustration, this

y[4][3] /. ParametricNDSolve[{
y''[x] + y'[x] + y[x] == 0,
y'[10] == y[10],
y[0] == u
}, y, {x, 0, 10}, u]


will evaluate just fine, but changing the y' initial condition to

y[4][3] /. ParametricNDSolve[{
y''[x] + y'[x] + y[x] == 0,
y'[u] == y[u],
y[0] == u
}, y, {x, 0, 10}, u]


will return the error

ParametricNDSolve::ndsv: Cannot find starting value for the variable y'.


In your case, of course, you've got your boundary conditions and that's that, and you can't really change them without changing your problem. What you can do is either

• move away from ParametricNDSolve as a solver, in favour of functionalized constructs such as

f[u_?NumericQ] := f[u] = y /. First[NDSolve[{
y''[x] + y'[x] + y[x] == 0,
y'[u] == y[u],
y[0] == u
}, y, {x, 0, 10}]]


which evaluates f[4][3] just fine, or

• re-scale your problem to a new independent variable ζ with a fixed range, such as ζ=z/(z0-a-zb), for which this shouldn't be a problem.

• Thanks a lot for your answer! I have some very basic doubts. Why the "First" is so crucial here? Also if I have two such parameter how to plot the function as I did in my question? Nov 10, 2015 at 12:26
• The First is because NDSolve returns its results as a list of lists of replacement rules (as a parallel to DSolve, which can return multiple independent solutions, cf. the output of DSolve[y'[x] + x y'[x]^2 == 1, y, x]), as opposed to ParametricNDSolve, which returns a simple list of replacement rules. If you have more than one parameter simply use f[u_?NumericQ,v_?NumericQ] := f[u,v] = (whatever). For more on _?NumericQ, see here. Nov 10, 2015 at 13:21
• Got it now. Another quick confusion : Plot[Table[f[u][x],{u,1,10}],{x,0,10}] gives back all curves but with same color. Any way out how to get different colors for different curves? Nov 10, 2015 at 13:39
• Wrap the Table inside an Evaluate. This will force the Table to evaluate to an explicit list, which the Plot command then recognizes as things that need to be plotted in different styles. It's admittedly a bit of a bizarre idiom. Nov 10, 2015 at 14:52
• In my original question if I use the following code (following your suggestion), it doesn't give me the output. I must be making some silly mistake here! f[w_?NumericQ, z0_?NumericQ] := f[w, z0] = x /. First[NDSolve[{DE, x'[z0-a] == -I (z0 w)/(6 a^2) x[z0-a], x[zb] == 1},x ,{z,zb,z0-a}]] f[0.001, 2][zb] Nov 10, 2015 at 19:29