Infinite Expression Error from NDSolve

Question

I am trying to solve a differential equation numerically. So I have

zSolutionRule = NDSolve[
   {
    z''[x]*z[x] - z'[x]^2 - z[x]^2 - p0[x]*z[x]^3 == 0,
    p0'[x] == 0,
    z'[0]*E == z'[1],
    z[0] == 1,
    z[1] == E
   },
   {z, p0}, {x, 0, 1}
  ];

When I run this, though, I get an "infinite expression $\frac{1}{0}$ encountered" error, but I want to get the solution to the differential equation. How do I get NDSolve to give me the solution like it does with other differential equations? (I am using Mathematica 9.0.0.0.)

Background

At first I solved a related differential equation

ySolutionRule = NDSolve[
   {y''[x] - 1 - p0[x]*Exp[y[x]] == 0,
    p0'[x] == 0,
    y'[0] == y'[1],
    y[0] == 0,
    y[1] == 1
    },
   {y, p0}, {x, 0, 1}
   ];
ySolution = y /. ySolutionRule[[1, 1]];
chargeNormalization = (p0 /. ySolutionRule[[1, 2]])[.5];

and this worked just fine, but I figured I could try to change variables to see if eliminating the exponential would change how good the solution is if I change boundary conditions when y[1] is big.

So I rewrote the equation in terms of $z=e^y$. I can check that the exponential of ySolution satisfies the transformed differential equation and boundary conditions:

zSol[x_] := Exp[ySolution[x]];
zSol[1]/zSol[0] - E
zSol[0] - 1
zSol[1] - E

and

Plot[zSol''[x]*zSol[x] - zSol'[x]^2 - zSol[x]^2 -chargeNormalization*zSol[x]^3, {x, 0, 1}]

It does, up to some accuracy. So then I tell it to find my solution but instead, it gives me an infinite expression error.

Answer 1

By default, NDSolve uses a "shooting" method to satisfy the boundary conditions: it picks an initial condition and then evolves the equation to see what boundary values are produced by the initial conditions. It then rejiggers these initial conditions to produce the boundary values specified by the boundary conditions. Here is the Mathematica documentation for the "Shooting" Method. It turns out that, for my differential equation, NDSolve's initial guess for initial condition is a bad guess. The documentation says

The shooting method by default starts with zero initial conditions so that if there is a zero solution, it will be returned.

This should sound bad because $z=e^y$ so z[0]==0 does not seem like a good guess. So just tell NDSolve what initial guess to make:

zSolutionRule = NDSolve[
   {
    z''[x]*z[x] - z'[x]^2 - z[x]^2 - p0[x]*z[x]^3 == 0,
    p0'[x] == 0,
    z'[0]*E == z'[1],
    z[0] == 1,
    z[1] == E
    },
   {z, p0}, {x, 0, 1},
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {z[0] == 1, z'[0] == 1, 
   p0[0] == -.1}}
   ];
zSolution = z /. zSolutionRule[[1, 1]];
chargeNormalization = (p0 /. zSolutionRule[[1, 2]])[.5];

The above code does not give an error. We can verify that it gives a solution satisfying the boundary conditions:

zSolution[1]/zSolution[0] - E
zSolution[0] - 1
zSolution[1] - E

and

Plot[zSolution''[x]*zSolution[x] - zSolution'[x]^2 - zSolution[x]^2 - 
  chargeNormalization*zSolution[x]^3, {x, 0, 1}, PlotRange -> All]