Eigenvalue dependent boundary conditions- mathematica

I am dealing with an eigenvalue problem whose boundary conditions are also eigenvalue dependent.

Could anyone please comment whether Mathematica can numerically solve such a problem? For boundary condition independent of eigenvalues, I use NDEigenSystem.

A minimal working example is given here. The eigenvalue problem: $$-\frac{d^2 \psi}{dx^2} +x^2 \psi = E \psi$$ with two boundary conditions: $$\textrm{(i) }\psi = 0 \textrm{ at }x = 0$$ and $$\textrm{(ii) }\frac{d\psi}{dx}+E^2\psi = 0 \textrm{ at } x = 1$$ needs to be solved to calculate the eigenvalues, $$E$$, of this operator. This might seem to be a trivial task, but please be aware of the eigenvalue-dependent boundary condition. I would be very thankful if anybody could suggest how to solve such eigenvalue problem in mathematica.

• Yes it can..... – David G. Stork Apr 27 '19 at 4:01
• You won't get a certain answer without showing us the precise equations. – Henrik Schumacher Apr 27 '19 at 7:05
• There was a question where the eigenvalue was both in the ode and was where the bc was applied: mathematica.stackexchange.com/questions/166952/… – KraZug Apr 27 '19 at 13:15
• Have you tried using NDEigensystem? – KraZug Apr 28 '19 at 15:15

While we wait for anyone to use NDEigensystem to see whether that works for this case (I expect it does), I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions or the github for some more details.

First we install the package (only need to do this the first time):

Needs["PacletManager"]
PacletInstall["CompoundMatrixMethod",
"Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]

Then we first need to turn the ODEs into a matrix form $$\mathbf{y}'=\mathbf{A} \cdot \mathbf{y}$$, using my function ToMatrixSystem (note capital E is reserved, so use lower case here):

Needs["CompoundMatrixMethod"]

sys = ToMatrixSystem[-ψ''[x] + x^2 ψ[x] == e ψ[x],
{ψ[0] == 0, ψ'[1] + e^2 ψ[1] == 0}, ψ, {x, 0, 1}, e]

The object sys contains the matrix $$\mathbf{A}$$, as well as similar matrices for the boundary conditions and the range of integration.

Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $$\lambda$$; this is an analytic function whose roots coincide with eigenvalues of the original equation.

FindRoot will then find solutions for a given start point:

FindRoot[Evans[e, sys], {e, 1}]
(* {e -> 9.9609} *)

You can see there are a whole set of positive real eigenvalues:

Plot[Evans[e, sys], {e, 0, 500}]

There are also some imaginary roots (I think just the two):

FindRoot[Evans[e, sys], {e, I}]
(* {e -> 0.159709 + 1.09317 I} *)

ContourPlot[{Re[Evans[er + I ei , sys]] == 0, Im[Evans[er + I ei , sys]] == 0},
{er, -3, 3}, {ei, -3, 3}, PlotPoints -> 30]