How can I use FindRoot on an expression from NDSolve?

Question

I have a second order ODE that I can only solve numerically using NDSolve, but I then need to use the solution in FindRoot and am running into errors. A simplified but analogous problem is the following: Find the solution to $\phi^{\prime\prime}(u) = -\omega^2 \phi(u)$ on the interval $0<u<1$ with the boundary conditions $\phi(0)=\phi(1)=0$, which will only be true for certain values of $\omega$ ($\omega=n\pi$ in this case). The code I am using is

eqnp = p''[ u] + ω^2 p[u];

psol[ω_?NumericQ] = NDSolve[{eqnp == 0, p[0] == 0, p[1] == 0}, p, u];

FindRoot[psol[ω], {ω, 3}]

but I am receiving errors. How can I feed a solution from NDSolve into FindRoot? I know that in the specific case above everything can be done analytically, but my actual ODE must be solved numerically, so I need to figure out how to solve this simplified problem completely numerically.

EDIT: I have updated the post with the full problem. The ODE is both homogeneous and linear

λ = 0.00001; k = 1.0;
f[u_] := 1/(2λ u) (1-Sqrt[1-4λ+4λ u^2]);
A = Sqrt[1/2(1+Sqrt[1-4λ])];
K2[u_] := u^(-3/2)f[u](1 + 2λ u^2f'[u]);
K1[u_] := K2[u](4ω^2)/(A^2f[u]^2) - 4k^2u^(-1/2)(1 - λ(6u^2f'[u] + 4 u^3f''[u]));
eomu = 4u^3K2[u]ϕ''[u] + (6u^2K2[u] + 4u^3K2'[u])ϕ'[u] + K1[u]ϕ[u];

I define the boundary conditions at $u=1$ by finding the first few coefficients in the series solution (using normal Froebenius method). I then want to find the solution and use Dirichlet boundary conditions at $u=0$ to determine the eigenfrequencies. The first coefficient is $a1$ and the indicial exponent is $-\frac{I\omega}{2A}$ (the other solution to the indicial equation is unphysical and is discarded)

a1 = (k^2(1 + Sqrt[1 - 4λ])(1 + 8λ) + 2(-1 + 2λ)ω^2 + (I Sqrt[1 + Sqrt[1 - 4λ]] (12λ ω))/Sqrt[2])/(2Sqrt[2] Sqrt[1 + Sqrt[1 - 4λ]] (Sqrt[1 + Sqrt[1 - 4λ]]/Sqrt[2] - I ω))

(I have separate code to find an arbitrary number of coefficients, but to enter them here would be too much) Then define the solution to use as boundary conditions

ϕhnum[u_] := (1 - u)^(-I ω/(2A))(1 + a1(1 - u))
numh = 0.999999;
numb = 0.000001;

Then I use the code suggested by Jens

ϕsol[ω_?NumericQ] := ϕ /. First@NDSolve[{eomu == 0, ϕ[numh] == ϕhnum[numh], ϕ'[
    numh] == ϕhnum'[numh]}, ϕ, {u, numb, numh}];

FindRoot[ϕsol[ω][numb], {ω, 1.0 - 1.0 I}]

When I do this I get the following error

FindRoot::lstol: The line search decreased the step size to within tolerance specified     
by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the 
merit function. You may need more than MachinePrecision digits of working precision to 
meet these tolerances. >>

I know that for the first solution I expect something like $\omega = 1.95- 1.27 I$. Hopefully this clarifies the problem, and thanks again for any help.

Answer 1

There are several things that need to be corrected. First of all, to get a non-trivial solution from NDSolve you have to

1. specify the start and end of the interval for u
2. Specify two initial conditions for a second-order differential equation. Initial and final conditions (corresponding to your boundary-value problem) won't work.

Since you want Dirichlet boundary conditions, you can choose the second initial condition to be a non-vanishing derivative. With this, you'll get a solution that doesn't satisfy the Dirichlet boundary condition p[1] == 0.

Then the job of FindRoot will be to adjust $\omega$ until the latter condition is saitsfied:

eqnp = p''[u] + ω^2 p[u];

psol[ω_?NumericQ] := p /. First@NDSolve[{
      eqnp == 0,
      p[0] == 0,
      p'[0] == 1
      },
     p, {u, 0, 1}];

FindRoot[psol[ω][1], {ω, 3}]

(* ==> {ω -> 3.14159} *)

To make FindRoot work with the solution, I also had to convert the output of NDSolve to an actual function first. The output is a rule, of which I take the First part which then replaces the generic p. So in FindRoot, the argument psol[ω][1] is the solution function for a given $\omega$, evaluated at the right-hand boundary u=1. When that's zero, you've found the correct value of $\omega$.

Answer 2

Here's a more concise version of Jens's solution, using the new NDSolve-related features introduced in version 9:

fun = 
  ParametricNDSolveValue[{p''[u] + ω^2 p[u] == 0, p[0] == 0, p'[0] == 1}, 
      p[1], {u, 0, 1}, ω]

FindRoot[fun[ω], {ω, 1}]

(* ==> {ω -> 3.14159} *)

This will find ω so that p[1] == 0.

---

Taking the fixed equation from your last update, and applying the same method works:

parfun = ParametricNDSolveValue[{eomu == 
    0, \[Phi][numh] == \[Phi]hnum[numh], \[Phi]'[numh] == \[Phi]hnum'[
     numh]}, \[Phi], {u, numb, numh}, \[Omega]]

Attempting FindRoot gives the error you encounter and gives no correct solution. There's a good chance that there's no solution: there's no ω so that parfun[ω][numb] == 0