Number of ways of selling 5 products for 4 customers
This is a straightforward application of the Polya Enumeration Theorem. Suppose $Z(S_n)$ is the cycle index of the symmetric group given by the recurrence by Lovasz which is
$$Z(S_n) = \frac{1}{n} \sum_{l=1}^n a_l Z(S_{n-l}) \quad\text{where}\quad Z(S_0) = 1.$$
We obtain by PET the following formula for the OGF of multisets (of products) on $n$ slots (the customers)
$$Z(S_n)\left(\sum_{k=1}^m X_k\right).$$
Note however that we have a restriction on these multisets which is that there are at most $q_k$ products of type $k.$ Now the generating function
$$\prod_{k=1}^m \frac{1}{1-X_k} Z(S_n)\left(\sum_{k=1}^m X_k\right)$$
has the property that the coefficient on $X_k^q$ counts multisets with at most $q$ instances of product $k$ (as opposed to exactly $q$). Therefore the end result is (here the square bracket denotes coefficient extraction)
$$\left[\prod_{k=1}^m X_k^{q_k}\right] \prod_{k=1}^m \frac{1}{1-X_k} Z(S_n)\left(\sum_{k=1}^m X_k\right).$$
The following code implements this formula using Maple. The formula can be optimized for particular values of the quantities $q_k$ in various ways e.g. by
$$\text{replacing}\quad 1/(1-X_k) \quad\text{by}\quad 1+X_k+\cdots+X_k^{q_k}.$$
(In fact numeric experiments indicate that this is not always an improvement.)
pet_varinto_cind := proc(poly, ind) local subs1, subs2, polyvars, indvars, v, pot, res; res := ind; polyvars := indets(poly); indvars := indets(ind); for v in indvars do pot := op(1, v); subs1 := [seq(polyvars[k]=polyvars[k]^pot, k=1..nops(polyvars))]; subs2 := [v=subs(subs1, poly)]; res := subs(subs2, res); od; res; end; pet_cycleind_symm := proc(n) local p, s; option remember; if n=0 then return 1; fi; expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n)); end; V := proc(q, n) local gf, m, vp; m := nops(q); gf := mul(1/(1-X[k]), k=1..m)* pet_varinto_cind(add(X[k], k=1..m), pet_cycleind_symm(n)); for vp to m do gf := coeftayl(gf, X[vp]=0, q[vp]); od; gf; end; W := n -> V([seq(n, q=1..n)], n);
A session with this code looks like this:
> V([3,3,3,3], 5); 40 > V([3,3,3,4], 5); 43 > V([5,5,5,5], 5); 56
and we have confirmation of the values given by the OP.
Additional verification may be obtained by treating the special case of $n$ products with max capacity $q_k = n$ and $n$ customers. Here we do not need the factor in front that we added because the cycle index itself enforces the limit of $n$ products of each type. We thus get the formula (all of whose different multinomials contribute one to the total):
$$\left. Z(S_n)\left(\sum_{k=1}^m X_k\right)\right|_{X_k=1}$$
Now recall the OGF of the multiset operator $\mathfrak{M}_{=n}$ which is $$Z(S_n) = [z^n] \exp\left(a_1 z + a_2 \frac{z^2}{2} + a_3 \frac{z^3}{3} + a_4 \frac{z^4}{4} +\cdots \right).$$
On substituting into the cycle index we put $$a_p = \left. \sum_{k=1}^m X_k^p\right|_{X_k=1} = n$$ because we have $m=n$ in this special case.
This yields the formula $$[z^n] \exp\left(n z + n \frac{z^2}{2} + n \frac{z^3}{3} + n \frac{z^4}{4} +\cdots \right) = [z^n] \exp\left(n\log\frac{1}{1-z}\right) \\ = [z^n] \frac{1}{(1-z)^n} = {2n-1\choose n-1}.$$
And indeed when we run the Maple code we obtain
> seq(W(n), n=1..18); 1, 3, 10, 35, 126, 462, 1716, 6435, 24310, 92378, 352716, 1352078, 5200300, 20058300, 77558760, 300540195, 1166803110, 4537567650 > seq(binomial(2*n-1,n-1), n=1..18); 1, 3, 10, 35, 126, 462, 1716, 6435, 24310, 92378, 352716, 1352078, 5200300, 20058300, 77558760, 300540195, 1166803110, 4537567650
This sequence is OEIS A001700.
The same computation works when all quantities $q_k$ of the $m$ products are at least $n$ where it yields the binomial coefficient $${n+m-1\choose m-1}$$ which appears to be an application of stars and bars.
Remark. The main formula above correctly represents the combinatorial givens of this problem but unlike in the example just presented with $n$ products, $n$ customers and capacity $n$ said formula computes all possible multisets and hence does not yield a reduction in complexity for the problem.
Addendum. We don't need the Polya Enumeration Theorem here. I will leave the above for the curious to consult. The formula is of course
$$[z^n] \prod_{k=1}^m (1+z+z^2+\cdots+z^{q_k}).$$
This yields for a uniform capacity $q$
$$[z^n] \prod_{k=1}^m \frac{1-z^{q+1}}{1-z} = [z^n] \frac{(1-z^{q+1})^m}{(1-z)^m} \\ = \sum_{p=0}^{\lfloor n/(q+1)\rfloor} {m\choose p} (-1)^p {m-1+n-p(q+1)\choose m-1}.$$
When $q\ge n$ only the term for $p=0$ contributes and we get the same formula as from the PET computation.