First, any strictly increasing linear function is a solution. Indeed, for any integer $k \geqslant 1$, the integer $a+k b$ divides $a^{2}+b \times(k a)=a(a+k b)$. Conversely, let's show that any solution is a linear function (which will be strictly increasing, since it takes values in $\mathbb{N}^{*}$).
Consider an integer $n \geqslant C$, and let $\varphi=f(1)$. Then

$$\varphi^{2}+f(n) \equiv n^{2}+f(n) \equiv 0 \quad(\bmod n+\varphi)$$

so $f(n)+\varphi^{2}$ is a non-zero multiple of $n+\varphi$, which means that $f(n) \geqslant n+\varphi-\varphi^{2}$.
On the other hand, since $1+f(n)$ divides $1+\varphi n$, let $g(n)=(1+\varphi n) /(1+f(n))$. Then

$$\varphi(1+f(n)) \geqslant \varphi\left(n+1+\varphi-\varphi^{2}\right)=(1+f(n)) g(n)-1+\left(1+\varphi-\varphi^{2}\right) \varphi$$

This means that $\Phi \geqslant(1+f(n))(g(n)-\varphi)$, where we have set $\Phi=\varphi^{3}-\varphi^{2}-\varphi+1$. Consequently, if $n \geqslant \Phi+\varphi^{2}-\varphi$, we know that $1+f(n)>\Phi$, and thus that $g(n) \leqslant \varphi$.
Now choose an integer $a \geqslant 1$, and let's prove that $f(a)=\varphi a$. To do this, we construct an integer $n \geqslant \max \left\{\Phi+\varphi^{2}-\varphi, C\right\}$ as follows:
$\triangleright$ we take $n \equiv 2(\bmod 4)$ if $\varphi$ is odd, and $n$ odd if $\varphi$ is even: in all cases, $n \varphi+1$ is odd;
$\triangleright$ for any prime $p \leqslant \max \{a, \varphi\}$ odd, we choose $n \equiv 1(\bmod p)$ or $n \equiv 2$ $(\bmod p)$ so that $n \varphi \not \equiv-1(\bmod p)$.
By the Chinese remainder theorem, such a choice is feasible, and for infinitely many integers $n$. Without loss of generality, we even assume that $n \geqslant 2 \max \{\varphi a, f(a)\}$.
Then $g(n)$ is a divisor of $1+\varphi n$, and we know that $g(n) \leqslant \varphi$. By construction, the integer $1+\varphi n$ has no prime factor $p \leqslant \varphi$, so $g(n)=1$, and thus $f(n)=\varphi n$. But then $a+\varphi n=a+f(n)$ divides

$$\left(a^{2}+n f(a)\right)+(a+\varphi n)(\varphi n-a)=(\varphi n)^{2}+n f(a)=\left(\varphi^{2} n+f(a)\right) n$$

However, we have constructed $n$ so that it has no common odd prime factor with $a$, and is not divisible by 4. Therefore, if we set $d=\operatorname{PGCD}(a+\varphi n, n)$, then $d=\operatorname{PGCD}(a, n)$ divides 2. Then, if we set $\alpha=(a+\varphi n) / d$ and $n^{\prime}=n / d$, we observe that $\operatorname{PGCD}\left(\alpha, n^{\prime}\right)=1$ and that $\alpha$ divides $\left(\varphi^{2} n+f(a)\right) n^{\prime}$, so $\alpha$ divides $\varphi^{2} n+f(a)$. The integer $\alpha$ therefore also divides

$$\left(\varphi^{2} n+f(a)\right)-d \varphi \alpha=f(a)-\varphi a$$

However, since $n \geqslant 2 \max \{\varphi a, f(a)\}$, we know that

$$\alpha>\varphi n / d \geqslant \varphi n / 2 \geqslant \max \{\varphi a, f(a)\} \geqslant|f(a)-\varphi a| .$$

Therefore, $f(a)=\varphi a$, which concludes.