Calculation of Bessel Functions

I want to calculate the Bessel function, given by

$$J_\alpha (\beta) = \sum_{m=0}^{\infty}\frac{(-1)^m}{m!\Gamma(m+\alpha +1)} \left(\frac{\beta}{2}\right)^{2m}$$

I know there are some tables that exist for this, but I want to keep the $\beta$ variable (i.e. I want a symbolic form in terms of $\beta$). If there is a way to simplify the summation part of the equation and leave an equation only in terms of $\beta$, that would be very helpful. (I see there is a dependence on $2m$, but I would like to see a way to break down the "other half" of the equation.)

Another question I have is: how is this calculated for $\beta$ values that are greater than $1$? It seems to me that this would give an infinite sum.

I am looking for something for $\alpha=1,3,5$ and $\beta=4$.

Thanks in advance.

As I alluded to in the comments, in general one would have to write a book chapter's worth of paragraphs to talk about the evaluation of Bessel functions for various argument ranges. Here, things are easier, since I only have to deal with integer orders of modest size. I shall now demonstrate one of my favorite methods, due to Yudell Luke.

Our starting point here is the pair of integrals

$$J_n(x)=\begin{cases}\frac2{\pi}\int_0^{\pi/2}\cos(x\cos\,u)\cos\,nu\;\mathrm du&n\text{ even}\\\frac2{\pi}\int_0^{\pi/2}\sin(x\sin\,u)\sin\,nu\;\mathrm du&n\text{ odd}\end{cases}$$

Two very useful methods for numerically evaluating these integrals are the trapezoidal rule and the midpoint rule. In a sense, these two are very accurate methods for the job, thanks to the Euler-Maclaurin formula. (See this for a deeper discussion.)

Using the odd order case as a concrete example, there is the following approximation which uses the (sadly lesser-known) midpoint rule:

$$J_n(x)\approx\frac1{m}\sum_{k=0}^{m-1}\sin\left(x\sin\left(\frac{\pi}{2m}\left(k+\frac12\right)\right)\right)\sin\left(\frac{\pi n}{2m}\left(k+\frac12\right)\right)$$

where $m$ is an appropriately chosen integer. For the particular case described in your question, taking $m=8$ gives approximations good to at least ten digits. Increase $m$ as needed.

In the case of even $n$, just replace all sines with cosines.

Again, this method is only suitable for modest integer values of $n$ and modest values of $x$; other methods might be more accurate, more efficient, or both for other argument ranges.

Here's another method you might want to consider, if you're in the business of generating a sequence of Bessel functions of fixed argument and consecutive integer orders. To describe the algorithm, which is due to J.C.P. Miller, we first take for granted the inequality ($x$ here is assumed real)

$$|J_n(x)|\leq\frac1{n!}\left|\frac{x}{2}\right|^n$$

and the series

$$1=J_0(x)+2\sum_{k=1}^\infty J_{2k}(x)$$

as well as the recurrence relation

$$\frac{J_n(x)}{J_{n-1}(x)}=\frac{x}{2n-x\frac{J_{n+1}(x)}{J_n(x)}}$$

Miller's idea is to first use an estimate like the inequality I gave to reckon an integer $n^\ast$ such that $\frac{J_{n^\ast}(x)}{J_{n^\ast-1}(x)}$ is smaller than machine epsilon. Having done so, pick some arbitrary value as a starting point (essentially, $f\,J_{n^\ast}(x)$ with $f$ an unknown constant), apply the recurrence backwards an appropriate number of times, while accumulating an unnormalized sum ($f\,J_0(x)+f\,J_1(x)+\cdots$). Once you've ended your recurrence, you can use the sum to normalize the recurrence values you stored along the way, which yields the Bessel function values you need.

To be more concrete, I shall present a Mathematica implementation of Miller's algorithm (which should be easily translatable to your favorite computing environment). I have chosen $n^\ast=24$ here; using the inequality with $x=4$, we have $|J_{24}(4)|\leq\frac{(4/2)^24}{24!}\approx 2.7\times10^{-17}$

x = N[4, 20];
n = 24;
(*hl accumulates ratios of Bessels h; s is the unnormalized sum*)
h = 0; s = 0; hl = {};
Do[
  h = x/(2 k - x h); (*recurrence relation*)
  hl = {h, hl};
  s = h (s + 2 Boole[EvenQ[k]]); , (*i.e., add 2 if k is even, else 0*)
  {k, n - 1, 1, -1}];
hl = Flatten[{1/(1 + s), hl}]; (*numerator is the value of the series*)
Do[hl[[k]] *= hl[[k - 1]], {k, 2, Length[hl]}];
hl

After executing the snippet, hl holds approximations to $J_0(4),J_1(4),\dots,J_{23}(4)$. When I tested it out, the first nineteen values generated were good to at least ten digits. Adapt the algorithm as needed.