I am looking for numerical method to solve the equation $$ \square \phi = \frac{\partial^2 \phi}{\partial x^2} - \frac{\partial^2 \phi}{\partial t^2} = \lambda \phi^3 \, $$ for real $\phi = \phi(x,t)$ and real $\lambda > 0$. I am interested in solutions on $x \in (0,l)$ and $t \in (0,t_{max})$ and supply the following boundary conditions: $\phi(0,t) = \phi(l,t) = 0$ and $\phi(x,0) = f(x)$ and $\dot{\phi}(x,0) = g(x)$.
I don't know much about numerical solutions of PDEs, but I have tried the most simple method I could think of and it was not converging to the correct solution.
Method I was thinking as follows. Let the step in $x$ and $t$ be the same, $\Delta$. Then, for a particular $x,t$ I can look at the partial derivatives and write them as $$ \square \phi (x,t) \approx \frac{\phi(x+\Delta,t)+\phi(x-\Delta,t)-\phi(x,t+\Delta)-\phi(x,t-\Delta)}{\Delta^2} \, , $$ i.e. it does not depend on the value at the specific point, but only on its neighbours. So, I created a 2D grid, assigned a random number to each point and then went as follows: For each point, I calculated $\square \phi$ and then set $\phi(x,t) = \left ( \frac{\square \phi(x,t)}{\lambda} \right )^{1/3}$. I then repeated the process several times. Total error, measured as sum of squares of differences between left and right hand sides of the equation, was indeed decreasing, but the results were completely unphysical (as in widely oscillating).
Could you please help me out with this?
SSF