Matlab numerical solving of a Second order Ode

Question

On matlab I want to solve numerically and plot this ode on $x\in [0,1]$ $$ g''(x)h(x)+g(x)(1-g(x)) = 0 $$ With (assume the 2 conditions) $a,b\in \mathbb R,\ g(0.5)=a$ and $g'(0.5)=b$ and where $h$ is given by $$ h(x) = \frac{1}{2\pi}e^{{-(\phi^{-1}(x)})^2} $$ and $\phi^{-1}$ is the inverse gaussian cumulative distribution.

On matlab it exists the function norminv for $\phi^{-1}$ (doc).

So I have to solve backwardly and forwardly the ODE (with an euler method for instance).

So as a beginner I have mainly 2 problems :

1. I do not know how to use norminv (which is not a symbolic function) inside odevectorfield and ode45 etc
2. How to tell to matlab to solve the ode backwardly (from $x=0.5$ to $x=0$) and forwardly (from $x=0.5$ to $x=1$)? and plot the final curve ?
3. does my below code is a good way to solve numerically and plot this ODE ? or do I have to use other packages ? dsolve for instance ? If someone has a solution or a code or a closed example of what I want ?

my code (is just the idea of the algorithm) that does not work at all is :

syms x,h(x),y(x) 

h(x) := 1./(2.*pi)*exp(-norminv(x)^2)%problem is norminv is not symbolic function

[V] = odeToVectorField(diff(y(x), 2) == y(x)*(y(x)-1.)./h(x)) 

M = matlabFunction(V,'vars', {'x','Y'})

sol = ode45(M,[0 1],[0.5 b]); %on x in[0, 1] and condition is g'(0.5)=b but problem is backwardly and forwardly is not specified to matlab.

plot(sol); %wrong since I have not specified on [0, 1]

But it does not work at all.

---

Extra informations :

For $\epsilon>0$ small, I have backwardly on from $0.5$ to $0$: $$ g(x-\epsilon) = g(x) - \epsilon g'(x)\\ g'(x-\epsilon) = g'(x) - \epsilon h^{-1}(x)g(x)(g(x)-1)\\ $$ And forwardly, for $\delta>0$ small from $0.5$ to $1$, $$ g(x+\delta) = g(x) + \delta g'(x)\\ g'(x+\delta) = g'(x) + \delta h^{-1}(x)g(x)(g(x)-1)\\ $$

Answer 1

Something along these lines, as @Ian suggested?

function odenorminv
  options = odeset('Stats', 'on', 'maxstep', 0.001);
  g0 = [0.5, 0.5];
  sol = ode45(@f, [0.5, 0.99], g0, options);
  plot(sol.x, sol.y)
  grid on
  legend('g','g''','location','northwest')
end

function ydot = f(x,g)
  ydot = [g(2); -g(1)*(1-g(1))/h(x)];
end

function rv = h(x)
  rv = exp(-norminv(x).^2);
end

Answer 2

One can "unpack" this problem by changing the independent parameter as $s=\Phi^{-1}(x)$, so that $x=\Phi(s)$ and consequently $h(x)=\frac1{2\pi}e^{-s^2}=\phi(s)^2$ where $\Phi'=\phi$ is the density of the normal distribution.

Computing the derivatives of the reparametrized function $u(s)=g(\Phi(s))$ gives $u'(s)=g'(\Phi(s))\phi(s)$ and $$ u''(s)=g''(\Phi(s))\phi(s)^2+g'(\Phi(s))\phi'(s) =-u(s)(1-u(s))+u'(s)\frac{\phi'(s)}{\phi(s)}. $$ or $$ u''(s)+su'(s)+u(s)(1-u(s))=0. $$ This now is, in forward direction, a mechanical system with steadily increasing friction and a potential function $\frac12u^2-\frac13u^3$. The steady states resp. extrema of the potential function are a minimum at $u=0$ and a maximum at $u=1$.

Using the energy function $$E=\frac12u'^2+\frac12u^2-\frac13u^3,$$

• if the energy at time $s=0$ is smaller than the local maximum $\frac16$ of the potential energy at any position $u$ between $0$ and $1$, the system will fall down into the minimum.
• if the energy is substantially larger than $\frac16$, the solution will move across the maximum and fall down the potential, that is, diverge to $+\infty$.
• for some initial energy in-between, the solution will converge toward the unstable maximum. The solution rapidly has to get to values close to $1$, else the energy bleed-off due to the increasing friction will drive the energy below the level of the local maximum.

As the equation is symmetric about $s=0, u=0.5$ or $x=0.5, y=0.5$, the whole analysis can be extended to the first half of the integration interval by exchanging the role of $u=0$ and $u=1$. Symmetric continuation will give one solution connecting stable to unstable stationary points $u(-\infty)=0$ to $u(\infty)=1$, and many solutions connecting the unstable to stable stationary points $u(-\infty)=1$ to $u(\infty)=0$.