Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am looking to calculate the Bessel function of the first kind $J_o(\beta)$. I am using the formula (referenced from wikipedia) to accomplish this.
$$J_\alpha (\beta) = \sum_{m=0}^{\infty}\frac{(-1)^m}{m!\Gamma(m+\alpha +1)} \left(\frac{\beta}{2}\right)^{2m+\alpha}$$

I am also aware of that MATLAB has a function which can calculate the solution as well. The function call is $besselj(nu,Z)$. However, what I am finding is that there is a discrepancy between what my $for$ loop calculates and what MATLAB outputs.

Does anyone see why this might be happening? I have included my code for reference:

loop_var = 100;
beta = 0.120;
alpha = 0;

sum = 0;
amp = 0;

[matl,ierr]  =besselj(alpha,beta)

for m=0:loop_var

Thanks for all comments and suggestions.

EDIT: Thanks to Fabian and Ed for catching that error.

I suppose as a "followup" question, the wikipedia formula and that MATLAB function seem to only match for "small" values of $\beta$. After $\beta$ is great than approximately 10, there is some error between the two values. Does anyone then know how the Matlab number is obtained?

share|cite|improve this question
Shouldn't it be for m=0:loop_var (I'm not good Matlab programmer, its just a hunch). – Fabian Sep 4 '12 at 21:15
@Fabian Good catch. I've fixed that. – suzu Sep 4 '12 at 21:24
@suzu I've addressed your edited question in my answer. – Emily Sep 4 '12 at 21:29
@EdGorcenski Thanks to you as well! I've now edited my question a bit with a follow-up that I hope has a quick fix as well :) – suzu Sep 4 '12 at 21:32
Since you have a $\beta^m$ term in there (actually $\beta^{2m+\alpha}$, but we can ignore the factor of 2 and the $\alpha$), then what happens as $m$ becomes large is that $\beta^m$ blows up for $\beta > 1$. There is no simple fix for this. The easiest solution might be to solve the Bessel's differential equation with those parameters, using an appropriate ODE method. – Emily Sep 4 '12 at 21:42
up vote 4 down vote accepted
for m=loop_var

This code only executes the loop once.

for m=0:loop_var

This code is correct (notice the 0: after m=).

For your edited question...

From the MATLAB besselj documentation:


The besselj function uses a Fortran MEX-file to call a library developed by D.E. Amos [3] [4].

[3] Amos, D.E., "A Subroutine Package for Bessel Functions of a Complex Argument and Nonnegative Order," Sandia National Laboratory Report, SAND85-1018, May, 1985.

[4] Amos, D.E., "A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order," Trans. Math. Software, 1986.

The first paper can be found here.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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