PDF with infinite basis

I used my usual tricks to show that you can define a family of orthonormal functions via the following recurrence relation:$$f_{n}\left(x\right) = \left(-1\right)^{n-1} \cdot \sqrt{\frac{2}{n}} \cdot x \cdot f_{n-1}\left(x\right) + \sqrt{1-\frac{1}{n}} \cdot f_{n-2}\left(x\right) $$$$f_{0}\left(x\right) = \pi^{-\frac{1}{4}} \cdot e^{-\frac{x^2}{2}}, f_{1}\left(x\right) = - \pi^{-\frac{1}{4}} \cdot \sqrt{2} \cdot x \cdot e^{-\frac{x^2}{2}}$$This is great for probability density functions as $\int^{\infty}_{-\infty} \left(\sum c_k f_k\right)^2 = \sum c_k^2$, meaning that if the norm of the coefficient vector is 1, then the square of that combination is guaranteed to be a pdf.

The code is below, but I'd like to go through the solution method I used to derive this. First, I used the formula for the $n$-th moment of a gaussian to derive the inner products of every function of the form $x^n \cdot e^{-\frac{x^2}{2}}$. For appropriate sigma and mu, this is given by $\int^{\infty}_{-\infty} x^n \cdot e^{-\frac{x^2}{2}} \cdot x^n \cdot e^{-\frac{x^2}{2}} dx = \Gamma\left(\frac{i+j+1}{2}\right)\cdot \left(1-((i+j)|2)\right)$. I created an $N \times N$ matrix of these terms to get the inner-products they share with each other. I then applied the QR-reshaping code from the other post to solve:

u,s,v = np.linalg.svd(inner_products)
h = (u/np.sqrt(s))
Q,R = flipper(h,'tr',flipped=True)

$R$ is the set of desired desired coefficient vectors to create this basis. It is however ugly, and to avoid using my brain I did the following; first, I saw a square root 2 term, and I squared the matrix of coefficients. From here, I converted these representations to fractions. Seeing that this worked well, I then identified the integer term each row could be multiplied by to force them to be integers. With that removed, it was simply a matter of identifying patterns in the integer rows. Knowing that we wanted an efficient solution in the form of a recurrence relation, and seeing repeated values, I shifted those out and used the row before, this presented the solution that could be used in the code.

dt = 0.001
theta = 0.5
N = 1000
t = np.arange(-8,8,dt)
w = np.exp(-(t**2)*theta)/(np.sqrt(np.sqrt(np.pi)))
B = np.zeros(shape=(N,len(w)))
B[:2,:] = np.vstack((w,-np.sqrt(2)*t*w))
for k in range(2,N):
    B[k,:] = np.sqrt(1-(1/k))*B[k-2,:]
    if k%2:
        B[k,:] -= t*np.sqrt(2/k)*B[k-1,:]
    else:
        B[k,:] += t*np.sqrt(2/k)*B[k-1,:]
print(np.allclose(B @ B.T,np.eye(N)/dt))
f = np.sin(t)
b = B @ f * dt
plt.plot(t,B.T@b)