# Contour plot for ParametricNDSolve

I want to solve a PDE in $(x,t)$ as a function of two parameters, $h$ and $w$, and then plot the times at which the solution at some $x=x_0$ is equal to a set value, as a function of $h$ and $w$. To put it more explicitly, I have a PDE

$$F\left(u(x,t),u_{xx}(x,t),u_t(x,t),x,h,w\right)=0$$

subject to a couple boundary conditions.

I am trying to plot the locus of points $(h,w,t)$ such that

$$u(x_0,t;h,w)=u_0.$$

Here's what I tried, using ParametricNDSolve and ContourPlot3D:

sol = ParametricNDSolve[{D[u[x, t], t]
==(1 - h*UnitBox[(x - 6)/w])*u[x, t]*(1 - u[x, t]) + D[u[x, t], {x, 2}],
u[x, 0] == UnitBox[x],u[40, t] == u[-40, t] == 0}, u,
{x, -40, 40}, {t, 0, 20}, {h, w}];

ContourPlot3D[u[h, w][20, t] /. sol == 0.9, {h, 0.1, 0.9}, {w, 0.5, 10}, {t, 0, 20}]


Unfortunately, this results in a series of errors about "[some long expression involving $u$] is neither a list of replacement rules nor a valid dispatch table and so cannot be used for replacing."

What have I done wrong here?

You have another problem but this is more a calculation problem so i will not say much about it. Your problem can be solved by just substituting the function beforehand:

sol = u /.
ParametricNDSolve[{D[u[x, t],
t] == (1 - h*UnitBox[(x - 6)/w])*u[x, t]*(1 - u[x, t]) +
D[u[x, t], {x, 2}], u[x, 0] == UnitBox[x],
u[40, t] == u[-40, t] == 0}, u, {x, -40, 40}, {t, 0, 20}, {h, w}]


And you're done. You can call this simply by:

ContourPlot3D[sol[h, w][20, t], {h, 0.1, 0.9}, {w, 0.5, 10}, {t, 0, 20}]


But now you will notice, that mathematica want unbelievable many points to solve your PDE. I mean, we can see why. You have a really big Range of values. So i'd consider to descrease your precision/accuracy or something. Or Increase your MaxPoints to get the full power (which will need a loooong time to plot).