Denote $P(x, y)$ as the governing equation of the problem. From $P(x, 0)$, we get $f(xf(0)) = (f(0))^{2021}$. If $f(0) \neq 0$, as $x$ is arbitrary, $f$ must be a constant function, say $f(x) = c$. Then $c = yc + c^{2021}$ for $y \in \mathbb{R}$, $c = 0$, a contradiction. Hence, $f(0) = 0$.
If $f(a) = 0$ for some $a \neq 0$, then $P(x, a)$ leads to $f(x) = 0$ for all $x \in \mathbb{R}$. Clearly, this function satisfies the governing equation. In the following, assume $f(a) \neq 0$ whenever $a \neq 0$.
From $P(0, 1)$, we get $f(1) = (f(1))^{2021}$, and $f(1) = \pm 1$. If $f(1) = -1$, $P(1, 1)$ gives $f(0) = -2$, a contradiction. Thus, $f(1) = 1$.
Now $P(x, 1)$ gives
$$f(x+1) = f(x) + 1. \qquad \textcircled{1}$$
Compare $P(x, y)$ and $P(x+1, y)$, to obtain
$$\begin{aligned} & f(xf(y) + f(y) + y^{2021}) \stackrel{P(x+1,y)}{=} yf(x+1) \\ & \quad +(y)^{2021} \stackrel{\textcircled{1}}{=} yf(x) + y + f(y)^{2021} \stackrel{P(x,y)}{=} y + f(xf(y) + y^{2021}). \end{aligned}$$
As $f(y) \neq 0$, $xf(y) + y^{2021}$ can take any real value $z$, and thus
$$f(z + f(y)) = y + f(z). \qquad \textcircled{2}$$
In (2), letting $z = 0$, we get $f(f(y)) = y$. Then, letting $y = f(w)$ in (2), we arrive at the following equation
$$f(z + f(f(w))) = f(z + w) = f(z) + f(w)$$
for any $z, w \in \mathbb{R}$.
---
Note that $P(0, y)$ gives $f(y^{2021}) = (f(y))^{2021}$, and so $P(x, y)$ becomes
$$\begin{aligned} yf(x) + (f(y))^{2021} &= f(xf(y) + y^{2021}) \\ &= f(xf(y)) + f(y^{2021}) \\ &= f(xf(y)) + (f(y))^{2021}, \end{aligned}$$
which implies $yf(x) = f(xf(y))$. Let $y = f(z)$ and use $f(f(z)) = z$ to
derive
$$f(xz) = f(x)f(z).$$
In particular, $f$ maps positive reals to positive reals. By the well-known
properties of Cauchy's multiplicative functional equation, $f(x) = x$, which
clearly satisfies the governing equation.
We conclude that there are only two functions: $f(x) = 0$ and $f(x) = x$.