Demo 10.1.5

Recast the TPBVP as an IVP by specifying a two-dimensional state as the original state and the first derivative (detailed in next section)

function ode!(f, y, lambda, r)
  f[1] = y[2]
  f[2] = lambda/y[1]^2 - y[2]/r
end;

function bc!(g, y, lambda, r)
  g[1] = y(0)[2]
  g[2] = y(1)[1] - 1
end;

domain = (eps(), 1.0);
est = [1, 0];

bvp = BVProblem(ode!, bc!, est, domain, 0.6);
y = solve(bvp);

plot(y, label = [L"w", L"w'"], legend=:right,
     xlabel = L"r", ylabel = "solution",
     title = "Solution of MEMS problems for lambda = 0.6")