, 2020

Solution:
Since $P(0)=0$, we see that $P$ has no constant term. Let $Q(x)=\frac{P\left(x^{2}\right)-x}{x}$ be a polynomial with degree at most $4039$. From the given values of $P$, we see that $Q(k)=0$ and $Q(-k)=-2$ for $k=1,2,3, \ldots, 2020$.
Now, consider the polynomial $R(x)=Q(x+1)-Q(x)$, which has degree at most $4038$. Then $R$ has roots $-2020,-2019, \ldots,-2,1,2, \ldots, 2019$, so
$$R(x)=a(x+2020) \cdots(x+2)(x-1) \cdots(x-2019)$$
for some real number $a$. Using $R(0)+R(-1)=Q(1)-Q(-1)=2$ yields $a=-\frac{1}{2020!2019!}$, so
$$Q(2021)=R(2020)+Q(2020)=-\frac{4040 \cdots 2022 \cdot 2019 \cdots 1}{2020!2019!}+0=-\frac{1}{2021}\binom{4040}{2020}$$
It follows that $P\left(2021^{2}\right)=2021 Q(2021)+2021=2021-\binom{4040}{2020}$.

Solution:
By Lagrange interpolation,
$$P(x)=\sum_{k=0}^{2020} k \prod_{\substack{0 \leq j \leq 2020 \\ j \neq k}} \frac{\left(x-j^{2}\right)}{\left(k^{2}-j^{2}\right)}=\sum_{k=0}^{2020} \frac{2(-1)^{k} k}{(2020-k)!(2020+k)!(x-k^{2})} \prod_{j=0}^{2020}\left(x-j^{2}\right)$$
Therefore, by applying Pascal's identity multiple times, we get that
$$\begin{aligned} P\left(2021^{2}\right) & =\sum_{k=0}^{2020} \frac{4042!(-1)^{k} k}{(2021-k)!(2021+k)!} \\ & =\sum_{k=0}^{2020}\binom{4042}{2021-k}(-1)^{k} k \\ & =2021-\left(\sum_{k=0}^{2021}(-1)^{k+1} k\left(\binom{4041}{2021-k}+\binom{4041}{2020-k}\right)\right) \\ & =2021-\left(\sum_{k=1}^{2021}(-1)^{k+1}\binom{4041}{2021-k}\right) \\ & =2021-\left(\sum_{k=1}^{2021}(-1)^{k+1}\left(\binom{4040}{2021-k}+\binom{4040}{2020-k}\right)\right) \\ & =2021-\binom{4040}{2020} . \end{aligned}$$