Which approach is correct to find $a_n$'s for Legendre Equation solution with boundary conditions?

Imagine we have a solution to a differential equation of the form $$T(r,x)=\sum_{n=0}^\infty a_nr^nP_n(x)$$ where $P_n(x)$ is the legendre polynomial satisfying the legendre equation of index $2$.

We want to find $a_n$ when $T$ satisfies boundary condition $$T(1,x)=\frac{1}{2}(3x^2-1) \qquad -1 \leq x \leq 1$$

APPROACH ONE

For our function, call it in general $f(x)=T(1,r)$ we have $$f(x)=\sum_{n=0}^\infty a_n P_n(x)$$ Now multiply both sides by $P_m(x)$ and integrate from $-1$ to $1$ to have $$\int_{-1}^1f(x)P_m(x)dx=\sum_{n=0}^\infty a_n\int_{-1}^1P_n(x)P_m(x)dx$$ by orthogonality and changing the name of the index this gives $$\int_{-1}^1f(x)P_n(x)dx=\frac{2a_n}{2n+1}$$ so finally $$a_n =\frac{2n+1}{2}\int_{-1}^1f(x)P_n(x)dx$$ now we notice that $f(x)=\frac{1}{2}(3x^2-1)=P_2(x)$ and hence by orthogonality $$a_n=\frac{2n+1}{2}\frac{2\cdot 2+1}{2}=\frac{5(2n+1)}{4}$$

APPROACH TWO

Similar to the above $$\int_{-1}^1\sum_{n=0}^\infty a_nP_n(x)P_m(x)dx=\int_{-1}^1P_2(x)P_m(x)dx$$ so that $$a_n = \frac{\int_{-1}^1P_2(x)P_m(x)dx}{\int_{-1}^1P^2_m(x)dx}= \begin{cases} \hfill 0 \hfill & n\neq m \\ \hfill 1 \hfill & n=m \\ \end{cases}$$

which one is correct?

In your first approach, the implication of $a_n=\cdots$ is only valid for $n=2$. Otherwise, you've got different Legendre polynomials so by orthogonality $a_n=0$. Additionally, you have $\int_{-1}^1 P_2(x)^2\,dx=\frac{2\cdot 2+1}{2}$, but the correct result is the reciprocal of this. Hence $a_2=\frac{5}{2}\cdot\frac{2}{5}=1$ and $a_n=0$ otherwise, in agreement with approach 2.
In my first approach I was very close actually, substituting before $f(x)=P_2(x)$ we have $$a_n = \frac{2n+1}{2}\frac{2m+1}{2}\delta_{mn}$$ which agrees with the second approach as shown by Semiclassical