Enumerating number of solutions to an equation

First you calculate the number of solutions in non-negative integers without worrying about the upper bound of $10$ on each variable. This is a standard stars-and-bars problem, reasonably well explained in the Wikipedia article. Then you use the inclusion-exclusion principle to get rid of the unwanted solutions.

In this case the first step gives you a preliminary figure of $$\binom{32+4-1}{4-1}=\binom{35}3=6545\;.$$

Now count the number of solutions that make $x_1$ too big. This means that $x_1\ge 11$, so the excess over $11$ in $x_1$ plus the values of $x_2,x_3$, and $x_4$ must add up to $32-11=21$. Each of these unwanted solutions therefore corresponds to a solution in non-negative integers to the equation $y_1+y_2+y_3+y_4=21$, and there are $$\binom{21+4-1}{4-1}=\binom{24}3=2024$$ of those. In fact there are $2024$ unwanted solutions for each of the four variables, so our next approximation is $6545-4\cdot2024=-1551$ solutions.

Of course this obviously isn’t right. The problem is that some solutions exceed the cap of $10$ on more than one variable. Every solution that exceeds the cap on two variables was removed twice when we subtracted $4\cdot 2024$ and therefore has to be added back in. Consider a solution that has both $x_1$ and $x_2$ greater than $10$. Then the excess in $x_1$, the excess in $x_2$ and the values of $x_3$ and $x_4$ must sum to $32-2\cdot11=10$, so we’re essentially counting solutions in non-negative integers to the equation $y_1+y_2+y_3+y_4=10$, of which there are $$\binom{10+4-1}{4-1}=\binom{13}3=286\;.$$ There are $\binom42=6$ pairs of variables, so we must add back in $6\cdot286=1716$ to get a better approximation of $-1551+1716=165$ solutions.

It’s impossible for more than two variables to exceed their quotas, since $3\cdot11=33>32$. Thus, no further corrections are needed, and the final answer is $165$ solutions meeting the original boundary conditions.

If you consider the following function $$ f_{\rm dim}(\epsilon)=\left(\frac{1-\epsilon^{11}}{1-\epsilon}\right)^{4}, $$ and expand at $\epsilon=0$ then coefficient in front of $\epsilon^{32}$ will give you the correct result, 165.

Explanation of why this works is given in my answer to This question.

The method can be obviously applied to generic case: Suppose there is equation $$\sum_{i=1}^n x_i=M$$ and we demand constraints $\lambda_i\leq x_i\leq \Lambda_i$. Question is how many solutions are there?. The answer is to consider

$$ f_{\rm dim}(\epsilon)=\prod_{i=1}^n\frac{\epsilon^{\lambda_i}-\epsilon^{\Lambda_i+1}}{1-\epsilon}\,, $$ expand this function at $\epsilon=0$ and find the coefficient of expansion at $\epsilon^{M}$.

Definitely, this method is a very efficient approach to use on a computer, much faster than generating all possible permutations as was suggested in another answer to your question.

For your particular example, this method can be also used to perform a computation by hands (though it might be not the case in generic situation). The required coefficient is given by the contour integral $\oint\frac{d\epsilon}{2\pi\,i}\frac{f_{\rm dim}(\epsilon)}{\epsilon^{33}}$ around the origin. But this integral can be also computed by residue at $\epsilon=\infty$ (note that at $\epsilon=1$ there is no pole). For the purpose of finding $1/\epsilon$ term in the large $\epsilon$ epxansion, the replacement $1-\epsilon^{11}\to-\epsilon^{11}$ can be used:

$$ \left(\frac{1-\epsilon^{11}}{1-\epsilon}\right)^{4}\frac 1{\epsilon^{33}}\simeq\epsilon^{11\times 4-33}\frac 1{(1-\epsilon)^4}=\epsilon^7\left(\frac 1{1-\epsilon^{-1}}\right)^4\to\epsilon^7\times\epsilon^{-8}\binom{8+4-1}{4-1}=\frac{165}{\epsilon} $$ After all, the answer is computed by a single Binomial coefficient $\binom{8+4-1}{4-1}$. This gives us a possibility to guess a nice trick. Consider a solution to the equation

$$ y_1+y_2+y_3+y_4=8 $$ with the only constraint $y_i\geq 0$. Then $x_i=10-y_i$ will be a solution to the original equation. It is easy to check that this is one to one map (with given boundary requirements), so $\binom{8+4-1}{4-1}$, the number of solutions for equation on $y$'s, is the desired answer.

In GAP, they can be computed via:

R:=RestrictedPartitions(32,[0..10],4);
S:=Union(List(R,r->Arrangements(r,4)));;
Size(S);

which gives 165.

Note that the first step generates unordered partitions of 32 into 4 parts, which I call $R$. Then I need to permute them in all possible ways, and take their union to create all possibilities, $S$.