# Worrisome result from pdetoode method in solving moving boundary condition problem

The problem I am trying to solve is the following :

$$(\frac{1}{c^2}\frac{\partial^2}{\partial t^2}-\frac{\partial^2}{\partial x^2}+1)\psi(x,t)=0$$

$$\psi (0,t)=\psi(L(t),t)=0$$

$$\psi(x,0)=\sqrt{\frac{1}{c\sqrt{1+\pi^2/4}}} sin(\pi \frac{x}{2})$$

$$\frac{d\psi (x,t)}{dt}\Bigr|_{t=0}=0$$

$$L(t)=2+sin(\omega t)$$

I have used the pdetoode method as was used here to solve this problem and here is my code (some codes might be redundant as I copied and modified the code from the link above) :

ydum = x/L[t];
eqn1 = 1/c^2 D[\[Psi][ydum, t], {t, 2}] -
D[\[Psi][ydum, t], {x, 2}] + (m^2 c^2)/h^2 \[Psi][ydum, t] ==
0 /. \[Psi][ydum, t] -> \[Psi][y, t] /. x -> y L[t] // Expand

m = 1;
c = 1;
h = 1;
\[Omega] = 1;
L[t_] := 2 + Sin[\[Omega] t];
bc = {\[Psi][0, t] == 0, \[Psi][1, t] == 0};
ic = {\[Psi][y, 0] ==
Sqrt[2 (m c)/(c Sqrt[m^2 c^2 + \[Pi]^2 h^2/L[0]^2] L[0])]
Sin[ \[Pi] y], D[\[Psi][y, t], t] == 0 /. t -> 0};
points = 25;
xdifforder = 4;
{yL, yR} = domain = {0, 1};
grid = Array[# &, points, domain];
pto = pdetoode[\[Psi][y, t], t, grid, xdifforder];
remove = #[[2 ;; -2]] &;
{ode, odebc, odeic} = MapAt[remove, pto /@ {eqn1, bc, ic}, 1];
soll = NDSolveValue[{ode, odebc, odeic[[2]], odeic[[1]]}, \[Psi] /@
grid, {t, 0, 50}];
sol1 = rebuild[#, grid, -1] &@soll
Manipulate[
Plot[sol1[y, t] /. y -> y L[t], {y, 0, 1},
PlotRange -> {-5, 10}], {t, 0, 50}]


The solution grows larger in time and this makes me anxious. Is this how it's supposed to be or is there something wrong with my code?

I have also solved this problem without the use of pdetoode (check out my previous question ) and I got different result :

I greatly appreciate anyone who would point out my error in both methods.

• The grid is too coarse, try e.g. points = 200; xdifforder = 2;. Also, it's better to observe a smaller time interval first e.g. {t, 0, 5}. Feb 19, 2022 at 13:36
• Thanks! Now the solution doesn't get larger in time. However, when I ramp up the $\omega$ to 1000, NDSolveValue spits out an error : NDSolveValue::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. How do I fix this? (By the way, what does xdifforder mean in pdetoode?) Feb 20, 2022 at 5:18
• 1. As to the icfail message, it's because you're currently discretize the system to a DAE system, to avoid it, you need to handle the b.c. in a more complicated way, see e.g. this post: mathematica.stackexchange.com/a/127411/1871 and this post mathematica.stackexchange.com/a/184285/1871 2. xdifforder sets the value of DifferenceOrder option of NDSolveFiniteDifferenceDerivative`, it's the order of difference formula for spatial discretization. Feb 20, 2022 at 14:31