Xrfgtjtk

Suppose you have an autonomous first order IVP $$x'(t)=f(x(t)), quad x(t_0)=x_0$$ where $x(t)inmathbb{R}$ and $x_0$ is a fixed point of the equation. Clearly, $$ x^0(t) equiv x_0$$ is a solution to this ODE. However, suppose there exists some other non-trivial solution $x^1(t)$ and that $f'(x(t))$ does not exist when $t=t_0$.

Is there a numerical procedure which will approximately solve the IVP and give $x^1(t)$ as a result?

Euler's method, $x(t+h) = x(t)+hf(x(t))%$, clearly converges to the trivial solution $x^0$, which can be seen via a simple induction argument. Other Taylor series methods are out since $f$ is not differentiable at $t_0$. It looks like basically any numerical procedure which converges to a solution of the the equation will converge to $x^0$.

Example: $x'(t) = sqrt{|x(t)|}, quad x(0)=0$

$x^0(t)equiv 0$ is one solution, and $x^1(t) = frac{1}{4}t^2text{sign}(t)$ is another (Grimshaw, "Nonlinear Ordinary Differential Equations", pg. 21 ex. 3). Any numerical procedure I throw at it converges to $x^0$, but I want one which approximates $x^1$.

Problems like this are to my knowledge often not discussed in numerical analysis. The problem you look at is clearly ill-posed.

However, there is the fairly new branch of mathematics called "probabilistic numerics" that tries to combine methods from statistics and numerical analysis. They in principal propose to reformulate the problem to a well-posed statistical one. This has been done for all kinds of numerical analysis, including the discretisation of ODEs. I suggest you look at Conrad et al. 2016 for probabilistic ODE solvers and maybe at http://www.probabilistic-numerics.org for the bigger image.

It looks like are trying explicit numerical methods. By their nature they can only give a single value at the new time -- multiple solutions are not obtained when sing explicit methods $y_{n+1} = g(y_n)$. Implicit methods involve solving an implicit equation for $y_{n+1}$ at each time-step. That should allow for multiple solutions since rootfinding methods must be able to handle multiple solutions.

Try Backward Euler (at first step): $y_1 - y_0 = dt; y_1^{1/3}$. Factoring on the left side yields $(1-dt;y_1^{-2/3})y_1=y_0$. We take $y_0=0$ and get two solutions for $y_1$: $y_1=0$ and $y_1=dt^{3/2}$. The second solution approximates the solution you want. Depending on the rootfinder you use, you might get either of these solutions.

In this case a semi-implicit solution that weights the right side evaluation as 1/3 at old time and 2/3 at new time will give zero numerical error on the first time-step for any step size - but that is specific to this example.