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.

  • 1
    $\begingroup$ Yes it can..... $\endgroup$ Apr 27, 2019 at 4:01
  • $\begingroup$ You won't get a certain answer without showing us the precise equations. $\endgroup$ Apr 27, 2019 at 7:05
  • $\begingroup$ There was a question where the eigenvalue was both in the ode and was where the bc was applied: mathematica.stackexchange.com/questions/166952/… $\endgroup$
    – SPPearce
    Apr 27, 2019 at 13:15
  • $\begingroup$ Have you tried using NDEigensystem? $\endgroup$
    – SPPearce
    Apr 28, 2019 at 15:15

1 Answer 1


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):

    "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):


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}]

enter image description here

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]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.