Problem with Legendre-Fourier series for sinx when the number of terms approaches infinity After I learned about Fourier series expansion, I understand orthogonality of trigonometric functions was the key when I calculate the coefficients of Fourier series. As I knew that Legendre Polynomials are also orthogonal with each other, I came up with a following problem:

Calculate Legendre-Fourier series for $f(x) = \sin{x}$


My solution
First, in this context, inner product of arbitrary functions $g, h$ is defined as $\langle g, h \rangle = \int_{-\pi}^{\pi} g(x)h(x) \mathrm{d}x$.
Next, define $Q_{n}(x)$ as an "adjusted" Legendre Polynomial of order $n$, which means $Q_{n}(x)$ satisfies $\langle Q_{m}(x), Q_{n}(x) \rangle = \delta_{mn}$. So $Q_{n}(x)$ can be written as
$$
Q_{n}(x) = \sqrt{\frac{2n+1}{2\pi}} P_{n}\left( \frac{x}{\pi} \right)
$$
where $P_{n}(x)$ is the Legendre Polynominal of order $n$.
Then $f(x)$ is expected to be expanded like
$$
f(x) = \sum_{n=0}^{\infty} c_n Q_{n}(x) \tag{1}\label{expected_expansion}
$$
where $c_{n}$ are the coefficients.
To calculate the coefficient $c_n$, I can utilize the orthogonality of $Q_{n}(x)$ by multiplying both sides of \eqref{expected_expansion} by $Q_{n}(x)$ and integrating them from $-\pi$ to $\pi$, and this derives
$$
c_n = \int_{-\pi}^{\pi} f(x) Q_{n}(x) \mathrm{d}x
$$
Thus, Legendre-Fourier series for $f(x) = \sin{x}$ can be expressed as
$$
\sum_{n=0}^{\infty} \left( \int_{-\pi}^{\pi} Q_{n}(x) \sin{x} \mathrm{d}x \right) Q_{n}(x) 
$$
$$\tag*{$\blacksquare$}$$
Issue
In order to confirm my answer, I defined $\widetilde{f}_{N}(x)$ as
$$
\widetilde{f}_{N}(x) = \sum_{n=0}^{N} \left( \int_{-\pi}^{\pi} Q_{n}(x) \sin{x} \mathrm{d}x \right) Q_{n}
$$
and plotted $\widetilde{f}_{N}(x)$ in Mathematica to see how well $\widetilde{f}_{N}(x)$ approximates $\sin{x}$ over $-\pi \le x \le \pi$ as $N$ increases.
However, although $\widetilde{f}_{N}(x)$ seems to fit $\sin{x}$ nicely when $N \le 20$, it suddenly fluctuates and diverges when $N \ge 21$. (Please look at this notebook) and I have no idea why this is happening.
Since $\sin{x}$ cannot be expressed as a finite sum of polynomials, I think $N$ should approach $\infty$ in order for $\widetilde{f}_{N}(x)$ to be equal to $f(x)$. Are there any errors in my solution or is this a kind of a bug of Mathematica?
Any help would be much appreciated. (Sorry for poor English)
Edit:
Thanks to the @Sergei Lytkin's comment, I now understand the computation of factorial of large numbers was causing this problem. However, even after I modified the notebook so that it uses recurrence formula to calculate Legendre Polynomial, the problem wouldn't vanish. Does anyone have other idea? Thanks in advance.
 A: In this answer, the coefficients are computed analytically.
We have $\sin\pi x=\sum_{n=1}^\infty (4n-1)a_n P_{2n-1}(x)$ with $a_n=\int_0^1 P_{2n-1}(x)\sin\pi x\,dx$.
Using $(2n+1)P_n(x)=\big[P_{n+1}(x)-P_{n-1}(x)\big]'$ and integration by parts, we get
$$\left.\begin{aligned}
C_n&:=\int_{-1}^1 P_n(x)\cos\pi x\,dx
\\S_n&:=\int_{-1}^1 P_n(x)\sin\pi x\,dx
\end{aligned}\right\}
\implies
\left\{\begin{aligned}
(2n+1)C_n&=\pi(S_{n+1}-S_{n-1})
\\(2n+1)S_n&=\pi(C_{n-1}-C_{n+1})
\end{aligned}\right.$$
which gives a computation of $a_n=S_{2n-1}/2$, with $C_0=C_1=S_0=0$ and $S_1=2/\pi$.
Alternatively, we have the following expression: $$a_n=\sum_{k=0}^{n-1}\frac{(-1)^k}{(2k)!}\frac{(2n+2k-1)!}{(2n-2k-1)!}\frac{1}{2^{2k}\pi^{2k+1}}.$$ A (long) way to get it: take $a_n=\frac12\int_{-1}^1$ and use Rodrigues' formula for $P_{2n-1}(x)$, then integrate by parts $2n-1$ times, then use Poisson's integral for Bessel functions to get $$a_n=\frac{(-1)^{n-1}}{\sqrt2}J_{2n-\frac12}(\pi)=(-1)^{n-1}j_{2n-1}(\pi),$$ and finally the explicit expression for the spherical Bessel function.
