we can find some other solutions such the zero function or the functions under the form $f(n)=a^2n$ or the constant functions under the form $f(n)=\dfrac{a^2}{2}$, or an other function $f$ which can be determined following this method.
Here I am just giving a way to build a function with this property
suppose that $a, b \in \mathbb{N}^{+} $ such that their sum is a perfect square, that is there exists $n \in \mathbb{N}_{≥2}$ such that
$$ a+b=n^2 $$
which is equivalent to
$$ a=n^2-b $$
which implies that
$$ f(a)+f(b)=f(n^2-b)+f(b) $$
then from here we can see that we are looking for a function $f$ such that
$$ \forall n \in \mathbb{N}_{\geq 2}, \forall b \in { 1, 2,..., [
\dfrac{n^2}{2}] }, \exists a(b,n) \in \mathbb{N}; f(n^2-b)+f(b)=(a(b,n))^2$$
Let's start building
for $n=2$ we need to give values to $f(b) $ and $a(b, 2)$ for $ b \in \{ 1, 2 \}$ such that
$$ f(4-b)+f(b)=(a(b,2))^2$$
that is we need to give values to $f(1), f(2), f(3)$ and $a(1, 2), a(2,2)$ such that
$$ \begin{cases} f(3)+f(1)=(a(1,2))^2 \\
f(2)+f(2)=(a(2,2))^2
\end{cases} $$
Then we find that we can associate a random values to $f(1), a(1,2), a(2,2)$ and we define $f(2), f(3)$ as follows
$$ \begin{cases} f(2)=\dfrac{(a(2,2))^2}{2} \\
f(3)=(a(1,2))^2-f(1) \end{cases}$$
before moving to the next step we need to note that if $n$ is even and by taking $ b= \left[ \dfrac{n^2}{2}\right]$ the equality $ f(n^2-b)+f(b)=(a(b,n))^2 $ gives
$$ f( [ \dfrac{n^2}{2}] )=\dfrac{a( [ \dfrac{n^2}{2}], n)^2}{2}...(*) $$
and we can choose any natural value for $a( [ \dfrac{n^2}{2}], n)$ to define $f( [ \dfrac{n^2}{2}] )$.
in the next steps we have to take in account that $f( [ \dfrac{n^2}{2}])$ has been already determined.
for $n=3$ we need to give values to $f(b) $ and $a(b, 2)$ for $ b \in \{ 1, 2, 3, 4\}$ such that
$$ f(9-b)+f(b)=(a(b,3))^2$$
that is
$$ \begin{cases} f(8)+f(1)=\dfrac{(a(1,3))^2}{2} \\ f(7)+f(2)=(a(2,3))^2\\ f(6)+f(3)=(a(3,3))^2\\ f(5)+f(4)=(a(4,3))^2
\end{cases}$$
Here we need to see that $f(2)$ and $f(8)$ have been already choosed from the equality $(*) $ that implies that $f(1)$ is not random as we find in the first step so we need to change the choice of $f(1)$ to be under the form
$$ f(1)=\dfrac{(a(1,3))^2}{2}-f(8)$$
and since $f(2), f(3)$ are well determined we can easily associate values to $f(7), f(6)$.
but $f(4), f(5)$ remain not determined, it's clear that we can associate
a random value to $f(4)$ and then we get $f(5)=(a(4,3))^2-f(4)$, but is this true?
Maybe in the next steps we will need to change one of these values as we have done with $f(1)$, and we iterate counting the next values of the function $f$.