#!/usr/bin/python3
# An unstable algorithm: spherical Bessel functions of the first kind from recursion

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import spherical_jn  # The j_n(x) function imported from
                                        # scipy, for comparison
                                        
def mybesselj(n, x):    # returns [j_0(x), j_1(x), ..., j_{n-1}(x)] computed by recursion
    values = np.zeros(n, dtype=float)  # table of function values
    values[0] = np.sin(x) / x                      # Given j_0(x) ...
    values[1] = np.sin(x) / x**2 - np.cos(x) / x   # ... and j_1(x) ...
    for k in range(2, n):                          # the others are computed recursively
        values[k] = (2*k - 1) / x * values[k-1] - values[k-2]
    return values

# Plot the results:
plt.plot(mybesselj(38, 10))  
for n in range(40):
    plt.plot(n, spherical_jn(n, 10), 'go')
plt.xlabel("n")
plt.ylabel("jn(10)")
plt.show()