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.

The next part of the program is to find the solution for the following equation.

\begin{equation} J_{0}(x)Y_{0}(x)=J_{2}(x)Y_{2}(x)\\ \end{equation}First of all, I constructed a new function \(f(x)=J_{2}(x)Y_{2}(x)-J_{0}(x)Y_{0}(x)\) and the question now becomes the root finding for \(f(x)=0\). The function \(f(x)\) is plotted in Fig. 3. Secondly, two root finding algorithms are constructed by two template routines called bisect and falsepo for bisection and false position methods(W. H. Press 2007) respectively. The templates take in an interval bracketed by \(x_{a}\) and \(x_{b}\), the calculation tolerance and the max number of iterations (say, \(10^{4}\) to prevent infinite loops). Given an interval, the program finds the root \(x_{r}\) within the interval to 7 significant digits such that the f