# solving Bessel function equation by hand

I have a Bessel function of the first kind given by the equation $$J_\alpha (\beta) = \sum_{m=0}^{\infty}\frac{(-1)^m}{m!\Gamma(m+\alpha +1)} \left(\frac{\beta}{2}\right)^{2m+\alpha}$$

I am trying to write a simple matlab code calculate the values of $J_{\alpha}(\beta)$. So far I have the following

beta = 4;
alpha = 1;
iteration = 3;

format long

for m = 1:iteration
J(m) = (((-1)^(m-1))/(factorial(m-1)*gamma(m+alpha)))*(beta/2)^(2*m+alpha-2);
end

out = sum(J)


For some reason, when I compare my results to the besselj() function in MS Excel, I do not get the correct answer. I have tried also making the number of summing terms higher (~30), but I'd like to find the fewest number of summing terms as possible. Is there something I am missing here in my code?

I suppose my other question is this. Intuitively, I thought there was an error in this formula, as I expected J to be a large number (in fact infinite) for $\beta$ values larger than 1. Why does this infinite summation work?

I am also aware there are other "forms" of this equation, in the form of integrals and some approximations, but I think this form is the simplest for me to work with right now.

Thanks in advance to all posters.

UPDATE: I used the provided documentation in Matlab, and tried using the formula of their implementation. However, I actually get the same answer as my formula (nothing exciting there), but still not the same answer as the besselj() function call. I'm stumped!

-
Any reason why you don't want to use Matlab's built in implementation of Bessel functions? (See mathworks.se/help/techdoc/ref/besselj.html) It is almost certainly both quicker and more accurate than what you could write yourself without working incredibly hard. – mrf Jun 7 '12 at 21:53
@mrf What you say is true and correct, but I'd like to write the whole summation as just a addition of a few terms if possible. eg from m = 0 to m = 3. Right now, I'm just trying to figure out how many terms would be considered enough to get a reasonable amount of accuracy, but I can't quite figure out what is wrong with the Matlab code that I have implemented above. To me, it seems like it should be possible to estimate the value of J with a finite number of terms, and hopefully a relatively small number of terms at that. – suzu Jun 7 '12 at 23:03
You forget to include the $m=0$ term in your sum. – mrf Jun 9 '12 at 20:59
@mrf, thanks. I just figured that out late last night and forgot to post here. It's the little things.... – suzu Jun 11 '12 at 4:42