Is it true that this function $f(n)=n^{13}$?

Let us denote by $f^{[k]}(n)$ the $k$th iterate of $f$. I cannot prove the claim, but I can prove that for all integers $n>1$ we have $$ \lim_{k\to\infty}\frac{\log f^{[k+1]}(n)}{\log f^{[k]}(n)}=13. $$ This is some kind of asymptotic evidence in favor of $f(n)=n^{13}$ being the only solution - alas, anything but conclusive.

This is seen as follows. We first prove that for all $k\ge3$ we have $$ f^{[k]}(n)=f^{[2]}(n)^{A_k} f(n)^{B_k} n^{C_k}\qquad(*) $$ for the sequence of vectors of positive integer determined by the recurrence relations $$ \left(\begin{array}{r} A_2\\B_2\\C_2\end{array}\right)=\left(\begin{array}{r} 1\\0\\0\end{array}\right),\qquad \left(\begin{array}{r} A_{k+1}\\B_{k+1}\\C_{k+1}\end{array}\right)=M\left(\begin{array}{r} A_k\\B_k\\C_k\end{array}\right), $$ where $M$ is the $3\times3$ matrix $$ M=\left(\begin{array}{crr} 1&1&0\\1&0&1\\2015&0&0\end{array}\right). $$ The proof follows from the given functional equation of $f$ by induction on $k$. The case $k=3$ is exactly the functional equation. The inductive step follows from the induction hypothesis by substituting $f(n)$ in place of $n$ and again applying the given functional equation.

The eigenvalues of $M$ are $\lambda_1=13$ and $\lambda_{2,3}=-6+i\sqrt{119}$. The key is that of these $\lambda_1$ has the largest absolute value. Furthermore, if we write the vector $$ (A_2,B_2,C_2)^T=x_1e_1+x_2e_2+x_3e_3 $$ in terms of unit eigenvectors $e_1,e_2,e_3$ belonging to the respective eigenvalues, we see that $x_1\neq0$.

For any $k\ge3$ we then have $$ (A_k,B_k,C_k)^T=\lambda_1^{k-2}x_1e_1+\lambda_2^{k-2}x_2e_2+\lambda_3^{k-2}x_3e_3. $$ For very large values of $k$ the first component dominates, and consequently $$ \lim_{k\to\infty}\frac{A_{k+1}}{A_k}=\lim_{k\to\infty}\frac{B_{k+1}}{B_k}=\lim_{k\to\infty}\frac{C_{k+1}}{C_k}=13. $$ It takes a while to see these limits if you calculate them. For some numerical support I fired up my Mathematica. The entrywise ratios of $M^{129}$ and $M^{128}$ are all in the interval $(12.9,13.2)$.

The claim follows from this by taking logarithms from $(*)$.

---

I don't know if this helps. It does seem to me that we should concentrate on large values of $n$ and asymptotics first. If only we could prove that $f$ must be a homomorphism of multiplicative monoids. Then it being strictly increasing would force it to a power function, and we know the exponent. If we know that $f$ is a power function the exponent can be determined without resorting to the above asymptotic gymnastics.

Re-arranging:

$h(n) = \dfrac{f(f(f(n)))}{f(f(n))\cdot f(n)} = n^{2015} $

Suppose $f(n) = n^{k}$:

$\dfrac{n^{k^3}}{n^{k^2 + k}} = n^{2015} $

$k^3 - k^2 - k - 2015 = 0$

which has solutions of $\{13, -6 \pm i \sqrt{119} \}$. The complex solutions oscillate, so $k = 13$. Clearly, $h(n)$ is unique and of the form $n^k$, and there's only one mapping from $f$ to $h$, so $f(n$) is unique.

EDIT: Regarding solutions not of the form $n^k$, define $g(n) = f(f(n))$

$f(g(n)) = g(f(n)) = g(n) \cdot f(n) \cdot n^{2015} $

We have to deal with the $n^{2015}$ term as it is of the form $n^k$. Suppose that $f(n) = \dfrac{l(n)}{n^{a}}$ and $g(n) = \dfrac{m(n)}{n^{b}}$ where $a+b = 2015$ and $l$ and $m$ are non-power series solutions by hypothesis:

$\dfrac{l\left(\dfrac{m(n)}{n^b}\right)}{n^a} = \dfrac{m(n)}{n^{b}} \cdot \dfrac{l(n)}{n^{a}} \cdot n^{2015}$

$l\left(\dfrac{m(n)}{n^b}\right) = m(n) \cdot l(n) \cdot n^{2015 - b}$

but $f(g(n)) = g(f(n))$, so

$\dfrac{m(l(n))}{n^b} = m(n) \cdot l(n) \cdot n^{2015 - b}$

$m(l(n)) = l(n) \cdot m(n) \cdot n^{2015} $

Which is what we started with. Therefore, both $f(g(n))$ and $g(f(n))$ must produce $n^{2015}$, but neither $f(n)$ nor $g(n)$ may contain $n^{\pm q}$ by hypothesis and then redundancy.

No, there are many such functions. For convenience (to avoid writing $\mathbb{N}^+$ too much) my intervals will only consist of natural numbers, so I'll write $[a, b]$ to mean $\mathbb{N}^+ \cap [a, b]$ and $(a, b)$ to mean $\mathbb{N}^+ \cap (a, b)$ below.

Consider any $f$ defined as follows: first let $f(1) = 1$ and $f(n) = n^{13}$ for $n$ of the form $n = 2^{13^k}$ ($k \geq 0$), so clearly $f$ satisfies the functional equation for $n = 1, 2^{13^k}$. We'll now inductively define $f$ on $A_k = [2^{13^k}, 2^{13^{k+1}}]$ for each $k \geq 0$, so that $f$ is strictly increasing on $A_k$. Note that since we already have $f(2^{13^k}) = 2^{13^{k+1}}$ and $f(2^{13^{k+1}}) = 2^{13^{k+2}}$, $f$ being strictly increasing will imply $f^{-1}(A_{k+1}) = A_k$.

First, for $k = 0, 1$, let $f|(2^{13^k}, 2^{13^{k+1}})$ be any strictly increasing function from $(2^{13^k}, 2^{13^{k+1}})$ to $(2^{13^{k+1}}, 2^{13^{k+2}})$ (where there is at least one such function since the second set is larger than the first), so clearly $f$ is increasing on $[2^{13^k}, 2^{13^{k+1}}] = A_k$.

Now let $k \geq 2$. Assume we've defined $f$ on $A_0, \dots, A_{k-1}$, hence on all $n \leq 2^{13^k}$. We define $f$ on $A_k$ in two steps. First, to ensure that $f$ satisfies the functional equation, we only need to correctly set the values of $f(f(f(n)))$, i.e. to correctly set the values of $f(a)$ for $a$ in the image of $f^2$. Since our $f$ will satisfy $f^{-1}(A_{i+1}) = A_i$ for each $i$, the image of $f^2$ in $A_k$ will be exactly $f(f(A_{k-2}))$; we already know these points since we've defined $f$ on $A_{k-2}$ and $A_{k-1}$. Second, since the definition in the first step will actually guarantee that $f$ is strictly increasing on the image of $f^2$, to make $f$ strictly increasing on all of $A_k$ we only need check that that there's enough room to assign $f$-values to the remaining points.

Step 1: Defining $f$ on $f(f(A_{k-2}))$. Write $f(f(A_{k-2})) = \{a_0, \dots, a_r\}$ where $a_0 < a_1 < \cdots < a_r$, so in particular since $f(f(2^{13^{k-2}})) = 2^{13^k}$ and $f(f(2^{13^{k-1}})) = 2^{13^{k+1}}$ we have $a_0 = 2^{13^k}$ and $a_r = 2^{13^{k+1}}$. For each $a_i$ we have $a_i = f(f(m_i))$ for some $m_i \in A_{k-2}$, hence for $0 < i < r$ we define $$f(a_i) := f(f(m_i)) f(m_i) m_i^{2015} = a_i f(m_i) m_i^{2015}$$ by necessity (to satisfy the functional equation), where this already holds for $i = 0, r$. Note that $m_0 < m_2 < \cdots < m_r$ since $f^2$ is strictly increasing on $A_{k-2}$.

Step 2: Defining $f$ on $A_k \setminus f(f(A_{k-2}))$. Having now defined all $f(a_i)$, we define $f$ on each $(a_i, a_{i+1})$ to be any strictly increasing function $(a_i, a_{i+1}) \to (f(a_i), f(a_{i+1}))$. There is more than one such function since \begin{align*} f(a_{i+1}) - f(a_i) &= a_{i+1} f(m_{i+1}) m_{i+1}^{2015} - a_i f(m_i) m_i^{2015} \\ &> (a_{i+1} - a_i) f(m_i) m_i^{2015} \\ &> a_{i+1} - a_i \end{align*} which means $|(f(a_i), f(a_{i+1}))| > |(a_i, a_{i+1})|$. $f$ is now defined and strictly increasing on each $[a_i, a_{i+1}]$, so it's defined and strictly increasing on $[a_0, a_r] = A_k$.

After the induction, we've defined an $f$ which is strictly increasing on each $A_k$, hence strictly increasing on all of $\mathbb{N}^+$. The functional equation holds, since for any $n \neq 1$, we have $n \in A_k$ for some $k$, hence by our definition of $f$ at $f(f(n)) \in A_{k+2}$, $f$ satisfies the functional equation at $n$. To conclude, note that for each $k$ we made at least one choice (since $|A_k| > |A_{k-2}|$, so the set $A_k \setminus f(f(A_{k-2}))$ is nonempty). This means there are uncountably many (at least $|2^{\mathbb{N}}|$) functions $f$ which can be constructed this way, each of which is strictly increasing and satisfies the functional equation.