Taylor made

We want the coefficient $C_N$ of $x^N$ in the series expansion (for $|x| < 1$ ) of $\begin{array}{rcl} f(x) & = & \displaystyle \frac{x^3}{(1-x)^3(1+x+x^2)} \; = \; \frac{x^3}{(1-x)^2(1-x^3)} \\ & = & \displaystyle \left(\sum_{n=0}^\infty (n+1)x^n\right)\left(\sum_{m=1}^\infty x^{3m}\right) \; = \; \sum_{n=0}^\infty \sum_{m=1}^\infty (n+1)x^{n+3m} \end{array}$ and hence $\begin{array}{rcl} C_N & = & \displaystyle \sum_{m=1}^{\lfloor \frac13N\rfloor} (N - 3m + 1) \\ & = & \displaystyle \lfloor \tfrac13N\rfloor (N+1) - \tfrac32\lfloor \tfrac13N\rfloor\big(\lfloor \tfrac13N\rfloor + 1\big) \\ & = & \displaystyle \tfrac12\lfloor \tfrac13N\rfloor\Big(2(N+1) - 3\big(\lfloor \tfrac13N\rfloor + 1\big)\Big) \\ & = & \displaystyle \tfrac12\lfloor \tfrac13N\rfloor\Big(2N - 1 - 3\lfloor \tfrac13N\rfloor\Big) \end{array}$ Since $2016 =3 \times672$ , the answer is $C_{2016} = \boxed{677040}$ .

I've used another approche:

Let $\omega$ and $\bar \omega$ be the complex roots of $1 + x + x^2 = 0$ $\begin{array}{rcl} 18f(x) & = & \displaystyle \frac{18x^3}{(1-x)^3(1+x+x^2)} \\ & = & \displaystyle -2(\omega -x)^{-1} -2(\bar \omega -x)^{-1} +4(1 -x)^{-1} -12(1 -x)^{-2} +6(1 -x)^{-3} \end{array}$

By recurrence we can see that: $\begin{array}{rcl} ({-2(\omega -x)^{-1}})^{(n)} & = & \displaystyle n!( -2(\omega -x)^{-(n+1)}) \\ ({-2(\bar \omega -x)^{-1}})^{(n)} & = & \displaystyle n!( -2(\bar \omega -x)^{-(n+1)}) \\ ({4(1 -x)^{-1}})^{(n)} & = & \displaystyle n!( 4(1 -x)^{-(n+1)}) \\ ({-12(1 -x)^{-2}})^{(n)} & = & \displaystyle n!( -12(n+1)(1 -x)^{-(n+2)}) \\ ({6(1 -x)^{-3}})^{(n)} & = & \displaystyle n!( 3(n+1)(n+2)(1 -x)^{-(n+3)}) \end{array}$

So that:

$\begin{array}{rcl} 18\frac {{f}^{(n)} (0)} {n!}& = & \displaystyle -2({\omega}^{-(n+1)} +{\bar \omega}^{-(n+1)}) + 4 -12(n+1) +3(n+1)(n+2) \\ & = & \displaystyle -2({\omega}^{-(n+1)} +{\bar \omega}^{-(n+1) } +1) + 6 -12(n+1) +3(n+1)(n+2) \\ & = & \displaystyle -2({\omega}^{-(n+1)} +{\bar \omega}^{-(n+1) } +1) + 3n(n-1) \end{array}$

knowing : $\begin{array}{rcl} {\omega}^{-(n+1)} +{\bar \omega}^{-(n+1) } + 1 & = & 3 \displaystyle \quad for \quad n \equiv 2 \pmod 3\\ & = & 0\quad otherwise \end{array}$

And finally: $\begin{array}{rcl} \displaystyle \frac {{f}^{(n)} (0)} {n!} & = & \displaystyle \frac{(n+1)(n-2)}{6} \quad for \quad n \equiv 2 \pmod 3 \\ & = & \displaystyle \frac{n(n-1)}{6} \quad otherwise \end{array}$

The solution I had in mind is similar to Mark's. For the sake of variety, let me write up a more informal version. Let's start with the geometric series $\frac{x^3}{1-x^3}=x^3+x^6+x^9+....$

If $A(x)=\sum_{n=0}^{\infty}a_nx^n$ then $\frac{A(x)}{1-x}=\sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}a_k\right)x^n$ : the coefficients of the power series of $\frac{A(x)}{1-x}$ are the partial sums of the coefficients of the power series of $A(x)$ . We will use this principle twice to get what we need: $\frac{x^3}{(1-x)(1-x^3)}=x^3+x^4+x^5+2x^6+2x^7+2x^8+3x^9+...+671x^{2015}+672x^{2016}+...$

(As Mark states, the coefficient of $x^n$ turns out to be the integer part of $\frac{n}{3}$ .)

Now the coefficient of $x^{2016}$ in $f(x)=\frac{x^3}{(1-x)^2(1-x^3)}=\frac{x^3}{(1-x)^3(1+x+x^2)}$ is $3\sum_{k=1}^{671}k+672=\boxed{677040}$

Let $g(x) = f(x).(1 + x + x^2)= ({\frac{x}{1-x}})^3$ $\hspace{3cm}$ then $g(x)=({\frac{1}{1-x} - 1})^3 = - \frac{1}{(x-1)^3} -1 -\frac{3}{(x-1)^2} -\frac{3}{x-1}$ $\hspace{10cm}$ Now differentiate $g(x)$ n times wrt $x$ and put $x=0$ , $\hspace{8cm}$ $g^n(0)=\frac{(n+2)!}{2}-3.(n+1)!+3.n!= \frac{(n-1)(n-2).n!}{2}$ $\hspace{3cm}$ and also, $g^n(x)=f^n(x)(1+x+x^2)+n(1+2x)f^{n-1}(x)+n(n-1)f^{n-2}(x)$ $\hspace{2cm}$ (see for yourselves!) $g^n(0)=f^n(0)+nf^{n-1}(0)+n(n-1)f^{n-2}(0)$ $\hspace{10cm}$ Further observe, $g^n(0)-ng^{n-1}(0)=f^n(0)-n(n-1)(n-2)f^{n-3}(0)=(n-2).n!$ $\hspace{4cm}$ Divide both sides by $n!$ to get, $\frac{f^n(0)}{n!}-\frac{f^{n-3}(0)}{(n-3)!}=n-2$
Sum up the above equation from $n=3$ to $n=2016$ to get
$\frac{f^{2016}(0)}{2016!}-\frac{f(0)}{0!}=1+4+7+...+2014$
$\frac{f^{2016}(0)}{2016!}=677040$