add comments on instability issue

src/CoulombTwoBody.jl

@@ -213,6 +213,9 @@ where ``\frac{1}{\mu} = \frac{1}{m_1}+\frac{1}{m_2}``, ``a_\mu = a_0 \frac{m_\ma
 !!! note
     The associated Laguerre polynomials $L_n^{k}(x)$, not the generalized Laguerre polynomials $L_n^{(\alpha)}(x)$, are used in this model.
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}
@@ -264,6 +267,9 @@ i^{|m|+m} \sqrt{\frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|} = (-1)^{\frac{|m|+m}{2}} \s
 @doc raw"""
 `legendre_polynomial(model::CoulombTwoBody, x; n=0, m=0)`
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}

src/HarmonicOscillator.jl

@@ -135,6 +135,9 @@ where ``\omega = \sqrt{k/m}``, ``\xi = \sqrt{\frac{m\omega}{\hbar}}x``, ``A_n =
 @doc raw"""
 `laguerre_polynomial(model::HarmonicOscillator, x; n=0)`
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}

src/HydrogenAtom.jl

@@ -179,6 +179,9 @@ where the Laguerre polynomials are defined as ``L_n(x) = \frac{1}{n!} \mathrm{e}
 !!! note
     The associated Laguerre polynomials $L_n^{k}(x)$, not the generalized Laguerre polynomials $L_n^{(\alpha)}(x)$, are used in this model.
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}
@@ -230,6 +233,9 @@ i^{|m|+m} \sqrt{\frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|} = (-1)^{\frac{|m|+m}{2}} \s
 @doc raw"""
 `legendre_polynomial(model::HydrogenAtom, x; n=0, m=0)`
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}

src/MorsePotential.jl

@@ -181,6 +181,9 @@ where ``\omega = \sqrt{k/µ}`` is defined.
 !!! note
     The generalized Laguerre polynomials $L_n^{(\alpha)}(x)$, not the associated Laguerre polynomials $L_n^{k}(x)$, are used in this model.
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}

src/PoschlTeller.jl

@@ -167,6 +167,9 @@ where ``\mu = \mu(n) = n_\mathrm{max}-n+1``, and ``n_\mathrm{max} = \left\lfloor
 @doc raw"""
 `legendre_polynomial(model::PoschlTeller, x; n=0, m=0)`
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Associated Legendre polynomials are the associated Legendre functions for integer indices. Here we use the same notation of the associated Legendre functions as in the model HydrogenAtom.
 
 ```math

src/RigidRotor.jl

@@ -150,6 +150,9 @@ i^{|m|+m} \sqrt{\frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|} = (-1)^{\frac{|m|+m}{2}} \s
 @doc raw"""
 `legendre_polynomial(model::RigidRotor, x; n=0, m=0)`
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}

src/SphericalOscillator.jl

@@ -174,6 +174,9 @@ where ``\gamma = \mu\omega/\hbar`` and ``\xi = \sqrt{\gamma}r = \sqrt{\mu\omega/
 !!! note
     The generalized Laguerre polynomials $L_n^{(\alpha)}(x)$, not the associated Laguerre polynomials $L_n^{k}(x)$, are used in this model.
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}
@@ -222,6 +225,9 @@ i^{|m|+m} \sqrt{\frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|} = (-1)^{\frac{|m|+m}{2}} \s
 @doc raw"""
 `legendre_polynomial(model::SphericalOscillator, x; n=0, m=0)`
 
+!!! note
+    The closed-form expression has been verified to be numerically stable up to $n=10$. For larger $n$, numerical instabilities may arise; a recurrence-relation implementation will address this limitation.
+
 Rodrigues' formula & closed-form:
 ```math
 \begin{aligned}