Assume strictly monotone increasing function; such that $f:N^{+}\to N^{+}$, $h$ for all $n\in N^{+}$, $$f(f(f(n)))=f(f(n))\cdot f(n)\cdot n^{2015}$$
Prove or disprove:$f(n)=n^{13}$
Put $n=1,f(1)=m$ $$f(f(m))=mf(m)$$ Put $n=m$, $$f(f(f(m)))=f(f(m))f(m)m^{2015}\Longrightarrow f(mf(m))=m^{2016}(f(m))^2$$ What about following?
On
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.
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.