$p$ is a prime througout the solution $P(1,1)\Longrightarrow 2(1)+1\mid 2f(1)+f(1)^2\mid 4f(1)+2f(1)^2$ but $2f(1)+1\mid f(1)(2f(1)+1)=2f(1)^2+f(1)$, so $2f(1)+1\mid 3f(1)\Longrightarrow f(1)=1$. $P(p,p)\Longrightarrow 2f(p)+p^2\mid 2p^2f(p)+f(p)^2\mid 4p^2f(p)+2f(p)^2$ but $2f(p)+p^2\mid f(p)(2f(p)+p^2)=2f(p)^2+p^2f(p)$, so $2f(p)+p^2\mid 3p^2f(p)\mid 6p^2f(p)$ but $2f(p)+p^2\mid 3p^2(2f(p)+p^2)$, so $2f(p)+p^2\mid 3p^4$. Now $2f(p)+p^2$ must be one of $3p^2, p^3, 3p^3, p^4, 3p^4$. $P(1,p)\Longrightarrow f(p)+p+1\mid p^2f(p)+f(p)+1$ but $f(p)+p+1\mid p^2(f(p)+p+1)=p^2f(p)+p^3+p^2$, so $f(p)+p+1\mid p^3+p^2-f(p)-1$, so $f(p)+p+1\mid p^3+p^2+p$. Now $2f(p)+p^2$ can only be $p^3$ or $3p^2$. Assume it's $p^3$, then $f(p)=\frac{p^3-p^2}{2}$. But $\frac{p^3-p^2}{2}+p+1\mid p^2+p^2+p\Longrightarrow \frac{p^3-p^2}{2}+p+1\mid 2p^2-p-2$, take huge $p$ give contradiction. So $f(p)=p^2$. Now $P(a,p)\Longrightarrow f(a)+ap+p^2\mid a^2f(a)+p^4+p^2f(a)$ but $f(a)+ap+p^2\mid (a^2-ap+p^2)(f(a)+ap+p^2)=a^2f(a)+p^4+p^2f(a)+a^3p-apf(a)$ so $f(a)+ap+p^2\mid a^3p-apf(a)$, take $p$ huge gives $f(a)=a^2$, which works.

why f(1)can't be 0

xyz123456 wrote: why f(1)can't be 0 $\mathbb{N}$ means positive integers in this problem.

person who has not done divisibility fe in two years forgets you can multiply divisibility statements before using euclidean algorithm The answer is $f(x)=x^2$ only, which works since $x^2+x+1 \mid x^4+x^2+1$. Let $P(a,b)$ denote the assertion Claim: For all odd primes $p$, we have $p^2 \mid f(p)$. Proof: Suppose $f(p)\neq p^2$, else we're already done. $P(p,p)$ yields $2f(p)+p^2 \mid f(p)(f(p)+2p^2) \implies 2f(p)+p^2 \mid f(p)(f(p)-p^2)$. If $p \nmid f(p)$, then $\gcd(2f(p),p^2,f(p))=1$, hence $2f(p)+p^2 \mid f(p)-p^2$, but this is a size contradiction. Thus $p \mid f(p)$ we can factor out a $p$ and write $2\tfrac{f(p)}{p}+p \mid f(p)(\tfrac{f(p)}{p}-p)$. Again, if $p \nmid \tfrac{f(p)}{p}$, then $\gcd(2\tfrac{f(p)}{p}+p,f(p))=\gcd(2\tfrac{f(p)}{p}+p,\tfrac{f(p)}{p})=1$, so $2\tfrac{f(p)}{p}+p \mid \tfrac{f(p)}{p}-p$ which is again a size contradiction. Thus $p^2 \mid f(p)$. Now from $P(1,1)$ we find $2f(1)+1 \mid 2f(1)+f(1)^2 \implies 2f(1)+1 \mid f(1)+2$, since $\gcd(2f(1)+1,f(1))=1$. For size reasons this means $f(1)=1$. Then, $P(p,1)$ (where $p$ is an odd prime) yields $$f(p)+p+1 \mid p^2f(p)+f(p)+1 \implies f(p)+p+1 \mid p(p^2+p+1).$$Since $p^2 \mid f(p)$, $p \nmid f(p)+p+1$, hence we require $f(p)+p+1 \mid p^2+p+1$. On the other hand, $f(p) \geq p^2$, so $f(p)=p^2$ for all primes $p$. Finally, let $a$ be arbitrary. For any odd prime $p$, we have $p^2+ap+f(a) \mid p^4+f(a)p^2+a^2f(a)$, which is equivalent to $p^2+ap+f(a) \mid p^2+\tfrac{f(a)}{a}p+f(a)$ by the Euclidean algorithm. Letting $p$ be large now implies $f(a)=a^2$, as desired. $\blacksquare$