I have been thinking about this problem for a few days already, out of curiosity. The number of non-negative integer solutions to $m_1x_1+\cdots+m_rx_r=n$ is equal to the $n$-th coefficient of
$$f(z)=\frac{1}{(1-x^{m_1})\cdots (1-x^{m_r})}$$
Right? So, it seems that partial fractions decomposition is the natural way to go. However, in general, the roots of $1-x^n$ are complex (roots of unity) and not very pretty. So, I'm sure that there is a natural link between this problem, cyclotomic polynomials and complex numbers that might give a better insight or be less tedious.
Can someone offer a systematic way to tackle problems like this? For a concrete example, how many non-negative integer solutions are there when solving $2x_1+3x_2+5x_3=2019$?
Edit: A related question on MSE I think something like the answer by Felix Martin is close to what I'm looking for. I'm working on it now to see if I can understand it. Meanwhile, any explanation or ideas are appreciated.
We can determine the number of solutions under these constraints using generating functions (indeed, using one of the main results of the study of generating functions):
If we interpret the equation using generating functions as: $$m_1x_1+\cdots+m_rx_r=n \quad\rightarrow\quad \sum_{x_1\ge 0}...\sum_{x_r\ge 0} x^{m_1x_1+...+m_rx_r}$$
Then the sum of all coefficients of $x^n$ with ${m_1x_1+...+m_rx_r}=n$ will be equal to the amount of solutions to the equation in the naturals.
We can do now start tansforming the generating function: $$ \sum_{x_1\ge 0}...\sum_{x_r\ge 0} x^{m_1x_1+...+m_rx_r} = \left(\sum_{x_1\ge0} x^{m_1 x_1}\right) ... \left(\sum_{x_r\ge0} x^{m_r x_r}\right) $$ Using the equivalence $$ \left(\sum_{x_i\ge0} x^{m_i x_i}\right) = \frac 1 {1-x^{m_i}} $$ we can get the generating function into a closed form: $$ \left(\sum_{x_1\ge0} x^{m_1 x_1}\right) ... \left(\sum_{x_r\ge0} x^{m_r x_r}\right) = \prod_{i=1}^r \frac 1 {1-x^{m_i}} $$
Using this we can get a closed formula for the sequence $(a_n)_{n\in\mathbb N}$, where $a_n$: $$ \sum_{x_1\ge 0}...\sum_{x_r\ge 0} x^{m_1x_1+...+m_rx_r}=\left(\sum_{x_1\ge0} x^{m_1 x_1}\right) ... \left(\sum_{x_r\ge0} x^{m_r x_r}\right) = \sum_{n\ge 0} a_n x^n $$ In other words, the closed formula will directly tell us the number of solutions to the initial equation in the naturals.
Now is, where the big result of generating functions comes in (actually, we only need a part of it):
If we have a generating function, and now it's form as fractal polynomial, i.e. if we have $$ \sum_{n\ge 0} a_n x^n = \frac{P(x)}{Q(x)}$$ where $P(x)$,$Q(x)$ are polynomials, then we can make the following statement of the closed formula of $a_n$:
Let $z_0,..,z_k$ be the reciprocal zeroes of $Q(x)$, and $v_0,..,v_k$ the multiplicity of these zeroes
(e.g. $Q(x)= x^2$ would have the zero $0$, with multiplicity $2$)
Then the closed form of $a_n$ has the form:
$$ a_n = \sum_{i=0}^k P_i(n) z_i^n $$ where $P_i$ is a polynomial of degree $<v_i$
Another way to find the solution is the following:
We know $$ \sum_{i\ge 0} a_i x^i = \prod_{i=1}^r \frac 1 {1-x^{m_i}} $$ , where $a_n$ is the number of solutions of the equation $m_1x_1+\cdots+m_rx_r=n $.
Then we can do the following: $$ \frac{d^n}{dx}\sum_{i\ge 0} a_i x^i = \sum_{i\ge n}(i)(i-1)...(i-n+1) a_i x^(i-n) $$ And therefore, if we substitute $x=0$ into the calculated differential, we directly get the number of solutions: $$ \frac{d^n}{dx}\sum_{i\ge 0} a_i x^i \overset{x=0}{=} a_n $$
Unless you can find a closed formula for the $n$-th differential of $\frac{P(x)}{Q(x)}$ this method however will get quite difficult for larger $n$.
Unluckily, these equations tend to get very complex very easily, so let's do it for an easy example (ironically, it's so easy you won't even need the method): $x_1 + 2x_2 = 200$
$$ x_1 + 2x_2 = 200 \rightarrow \sum_{x_1\ge 0}\sum_{x_2\ge 0} x^{x_1+2x_2} = \left(\sum_{x_1\ge 0} x^{x_1} \right) \left(\sum_{x_2\ge 0} x^{2x_2} \right) = \frac{1}{1-x}\frac{1}{1-x^2} = \frac 1 {x^3 - x^2 - x + 1} = \frac 1 {(x + 1)·(x - 1)^2} $$
Therefore, we have the zeroes $-1, 1$ with multiplicities $1,2$. And arrive at the closed formula $$ a_n = (1/(-1))^n (\alpha) + (1/(1))^n (\beta x + \gamma) = (-1)^n (\alpha) + (\beta n + \gamma) $$
We now manually calculate the first few solutions to our equation $x_1 + 2x_2 = n$ for $n=1,2,3$.
$|L(x_1 + 2x_2 = 1)| = 1 \quad (x_1 = 1) $
$|L(x_1 + 2x_2 = 2)| = 2 \quad (x_1 = 2 \lor x_2 = 1) $
$|L(x_1 + 2x_2 = 3)| = 2 \quad (x_1 = 3 \lor x_1 = 1, x_2=1)$
We arrive at the linear equation system: $$ \begin{align*} (-1)^0 (\alpha) + (\beta 0 + \gamma) = 1\\ (-1)^1 (\alpha) + (\beta 1 + \gamma) = 2\\ (-1)^2 (\alpha) + (\beta 2 + \gamma) = 2 \\ \end{align*} $$ Solving this and substituting the solution of $\alpha,\beta,\gamma$ into $a_n$ leads to: $$ a_n = (-1)^n·(1/4) + (1/2·n + 3/4) $$ Which is also the number of solutions to the equation $x_1+2x_2 = n$.
Therefore we have $|L(x_1+2x_2 = 200)| = a_200 = 101$.
Just as reminder that it for usual gets really difficult really quickly, here is the closed formula (before determining the coefficients) of your example:
(Here, the coefficients are the $c_i$)
And here it is, with coefficients solved: $$ a_n:=(√5/25 + 1/5)·COS(4·\pi·n/5) + (1/5 - √5/25)·COS(2·\pi·n/5) + COS(2·\pi·n/3)/9 - √3·SIN(2·\pi·n/3)/9 + COS(\pi·n)/8 + n^2/60 + n/6 + 131/360 - i·SIN(\pi·n)/8$$
And finally, therefore we have: $a_{2019} = 68276$, i.e. your initial linear equation has 68276 solutions in the naturals.