I am trying to solve the Eigenvalue problem $$ x^2 y''+ x y' + x^2 y = \lambda^2 y\,,\quad x\in(0,1)\,,\quad y(0)=0\,,\quad y'(1)=y(1) $$ The differential equation is the Bessel equation. The solution is given by the Bessel function of the 1st kind $y(x)=J_\lambda(x)$. Bessel functions of the 2nd kind is omitted because of their singularity at $x=0$ . The Eigenvalue $\lambda$ has to be determined using the boundary condition at $x=1$. Using the identity $$J'_\lambda(x)=\frac{1}{2}\left(J_{\lambda-1}\left(x\right)-J_{\lambda+1}\left(x\right)\right)$$ the boundary condition at $x=1$ can be rewritten as $$J_{\lambda-1}\left(1\right)-2J_{\lambda}\left(1\right)-J_{\lambda+1}\left(1\right)=0\,.$$ Now I would like to determine the $\lambda$ with the equation above using the fzero-routine of MATLAB/ GNU Octave. Code below.

%% initial guess
%% function handle
f=@(l) besselj(l-1,1)-2*besselj(l,1)+besselj(l+1,1)
%% set tolerance
%% determine zero

If you plot the function f, click here, you observe that most of the zeros are in $-\infty<x<0$.

To get a sorted sequence $\lambda_{k+1}<\lambda_k$, I would like to search from the "last" zero only in negative $x$-direction to the "next" one. Is there any way to implement this in the fzero-routine?

Thank in advance!

  • $\begingroup$ It's not clear for me what your question is. What do you mean with "from the "last" zero only in negative x-direction to the "next" one"? Normally, one just need a finite set of eigenvalues, then you can make a progressive search for this problem. $\endgroup$
    – nicoguaro
    Jan 13, 2015 at 19:51
  • $\begingroup$ Let me clarify the problem. I would like the first 10 eigenvalues. Let's say I found $\lambda_1$. Then $\lambda_2$ should have the property that there is no zero of the function f between $\lambda_1$ and $\lambda_2$. $\endgroup$ Jan 13, 2015 at 20:53
  • $\begingroup$ How can I program a progressive search effectively? How to give a step size? If I detect a change if the sign I could use bisection method, right? $\endgroup$ Jan 13, 2015 at 21:00
  • $\begingroup$ You can specify an interval in place of an initial guess for fzero to search in, e.g. xzero=fzero(f,[-1e-9,xzero_old],opts). $\endgroup$ Jan 13, 2015 at 22:16
  • $\begingroup$ Is there a sign error in your code in front of besselj(l+1,1)? $\endgroup$
    – Kirill
    Jan 13, 2015 at 22:35

2 Answers 2


As you can see in this plot of $\log|f(\lambda)|$, $$ f(\lambda) = J_{\lambda-1}(1) - 2J_\lambda(1) -J_{\lambda+1}(1), $$

enter image description here

the roots $\lambda_k$ are really regular, and are approximately equal to $-k$ (starting from $k\geq0$, $\lambda_0=1.23219$ is an exception).

So the way to get the $k$-th eigenvalue ($k\geq1$), is to bracket the root in $[-k-\frac12,-k+\frac12]$ and use fzero solver in that interval (it accepts a bracketing interval as the initial guess and stays within that interval).

Note that the magnitude of the function grows very quickly (something like $(2\lambda/e)^\lambda\sin\pi\lambda$), so once the function becomes quite large (for large $k$), it may be better to use $\lambda_k=-k$ as an approximation.

  • $\begingroup$ The code in my answer to a related question performs almost exactly this task. $\endgroup$ Jan 14, 2015 at 2:41
  • $\begingroup$ @DougLipinski Yes, but your approach requires you to guess a "comb" of locations to separate roots, and you might miss some roots easily. In this exact question it's possible to find all roots guaranteed. $\endgroup$
    – Kirill
    Jan 14, 2015 at 3:24
  • 1
    $\begingroup$ Indeed, just choose the intervals in my code to be the ones you suggested. I just wanted to point out some existing code that solves the OP's problem with the choice of intervals you suggested. You've shown the right "comb" for this problem. +1 from me. $\endgroup$ Jan 14, 2015 at 3:34
  • $\begingroup$ The bracket you proposed works quite well. Thanks a lot. The log-plot is quite helpful. $\endgroup$ Jan 14, 2015 at 10:15

Here is the code I used now.

%% function handle
f=@(l) besselj(l-1,1)-2*besselj(l,1)-besselj(l+1,1);
%% set tolerance
%% compute first positive eigenvalue
%% compute negative eigenvalues
for k=1:(kmax-1)
  if not(fl==1)
    error('iteration    did     not     convert.')

By the way, I would like to use this to verify my finite difference method.


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.