Solution (source code)

= Solution

Take each <partial derivative> while holding the other coordinate fixed. Differentiating the <polynomial> term by term gives
$$
\boxed{\begin{aligned}
f_x&=2x(y^4+y^2+1),\\
f_y&=4x^2y^3+3y^2+2x^2y-3,\\
f_{xx}&=2(y^4+y^2+1),\\
f_{xy}&=8xy^3+4xy,\\
f_{yy}&=12x^2y^2+6y+2x^2.
\end{aligned}}
$$
Since the <function> is a <polynomial>, both orders of the <mixed partial derivative> agree everywhere.