First Solution
It is the same as the coefficient of $3n$ in the expansion $(1+x+x^2+x^3+\cdots+x^{2n})^3$. The proof is left to the readers as an exercise.
Note that $(1+x+\cdots+x^{2n})^3 = \frac{(1-x^{2n+1})^3}{(1-x)^3} = (1 - 3x^{2n+1} + 3x^{4n+2} - x^{6n+3})(1-x)^{-3}$
But $(1-x)^{-3} = (1+3x + 6x^2 + 10x^3 + \cdots + \frac{(k+2)(k+1)}{2}x^k + cdots)$
So the coefficient of $3n$ in that expansion is:
$\frac{(3n+2)(3n+1)}{2} - 3\frac{n(n+1)}{2} = 3n^2 + 3n + 1$
Second Solution
We shall prove that there are $3n^2+3n+1 = (n+1)^3 - n^3$ ways, which has striking geometrical interpretations.
Create an equilateral triangle with altitude $3n$, and mark all the points that has integer distance to all the sides. These points will form a triangular lattice
It is a well-known fact that for any point in an equilateral triangle, the sum of distances to the sides is constant. Since the altitude of the triangle is $3n$, then for any point in that triangular lattice, the sum of distance to the sides is $3n$. Furthermore, the distances are integers. Thus, these points represent the number of ways to choose $3n$ balls from 3 colors.
However, the constraint that we only have $2n$ supply of each color is not yet enforced. We examine the points in the lattice whose distance to any of the sides is greater than $2n$, and we eliminate them. So we are left with a hexagonal lattice of side $n+1$. This is the projection of a half-shell of an $n+1$-cube that has $n$-cube carved out of it.
