The Bessel Functions of the first kind \(J_{\nu}(x)\) and the second kind \(Y_{\nu}(x)\) are two functions of interest in this program. The program incorporates two header files ”nr3.h” and ”bessel.h” and outputs the function values between \(x=0\) to \(x=20\). The first two Bessel functions \(J_{0}\), \(J_{1}\), \(Y_{0}\) and \(Y_{1}\) are built in within ”bessel.h” and we can use the recurrence formula to construct higher order Bessel functions.

\begin{equation} J_{\nu+1}(x)=\frac{2\nu}{x}J_{\nu}(x)-J_{\nu-1}(x)\\ \end{equation} \begin{equation} Y_{\nu+1}(x)=\frac{2\nu}{x}Y_{\nu}(x)-Y_{\nu-1}(x)\\ \end{equation}The first three \(J_{\nu}\) and \(Y_{\nu}\) are plotted by MATLAB as shown in Fig. 1 and Fig. 2.