There is no nice formula but there is a nice algorithm. The multiplication of polynomial generating functions mirrors the sequential drawing process. As such, computing with generating functions essentially performs an exhaustive enumeration of the possibilities. However, algebraic properties can be exploited to simplify the intermediate results (partial products), leading to a considerably streamlined calculation. These algebraic properties (formally: an ideal of polynomials) encode the two facts that (a) we are not concerned with how many of each ball is drawn, only which balls are drawn; and (b) it suffices to focus on the cases where not all balls are drawn.
This reply shows how these ideas can be used efficiently to compute an answer to the question--as well as to its obvious generalizations to more balls and different numbers and sizes of draws. It is followed by analyses, including a simulation, demonstrating the correctness of the answer. A final section discusses further issues that might be of interest, such as applying this approach to the simpler (unweighted) case.
To draw $k$ balls (with weighted probabilities) from an urn with $n$ balls we draw the first (but do not replace it), then the second (without replacement), and so on. Having observed the result, we replace all balls. Therefore the chance of seeing a sequence of balls with probabilities $p_1, p_2, \ldots, p_k$ equals
$$\Pr(p_1, p_2, \ldots, p_k) = p_1\frac{p_2}{1-p_1}\cdots\frac{p_k}{1-(p_1+p_2+\cdots+p_{k-1})}$$
and the chance of seeing a combination (unordered sequence) of balls with those probabilities is the sum over all permutations,
$$\Pr(\{p_1, p_2, \ldots, p_k\}) = \sum_{\sigma\in\mathcal{S}_k} \Pr(p_{\sigma(1)}, p_{\sigma(2)}, \ldots, p_{\sigma(k)}).$$
A simple and elegant way to handle these $\binom{n}{k}$ numbers is with a generating polynomial,
$$m_k(x_1, x_2, \ldots, x_n) = \sum_{i_1\lt i_2\lt \cdots \lt i_k}\Pr(\{p_{i_1}, p_{i_2}, \ldots, p_{i_k}\})x_{i_1}x_{i_2}\cdots x_{i_k}.$$
For instance, with probabilities $(4/10, 3/10, 1/10, 1/10, 1/10),$ we may compute
$$m_1(x_1, x_2, x_3, x_4, x_5) = \frac{4 x_1}{10}+\frac{3 x_2}{10}+\frac{x_3}{10}+\frac{x_4}{10}+\frac{x_5}{10},$$
$$m_2(x_1, x_2, x_3, x_4, x_5) = \frac{13 x_1 x_2}{35}+\frac{8 x_3 x_2}{105}+\frac{8 x_4 x_2}{105}+\cdots+\frac{x_3 x_5}{45}+\frac{x_4 x_5}{45},$$
and
$$m_3(x_1, x_2, x_3, x_4, x_5) =\frac{76}{315} x_1 x_2 x_3+\frac{76}{315} x_1 x_2 x_4\cdots+\frac{17}{504} x_2 x_4 x_5+\frac{1}{120} x_3 x_4 x_5.$$
Writing $x = (x_1, x_2, \ldots, x_n),$ the generating polynomial for the total outcome of a sequence of independent draws of sizes $k_1, k_2, \ldots, k_d$ is simply the product $m_{k_1}(x)m_{k_2}(x)\cdots m_{k_d}(x).$ Moreover, all we are interested in is whether all the $x_i$ appear in the outcome, not how many times each does. This allows the computation of the product to be simplified: at each multiplication we may throw away any monomial term that is a multiple of $x_1x_2\cdots x_n$ and replace any power of any $x_i$ by $x_i$ itself (which can be done by equating each $x_i^2$ with $x_i$ itself). The result will encode the chances of not obtaining all $n$ balls among these $d$ draws: the coefficient of each monomial $x_{i_1}x_{i_2}\cdots x_{i_s}$ will be the chance of obtaining exactly the balls $\{i_1, i_2, \ldots, i_s\}$. Setting all the $x_i=1$ will give the chance that not all balls are drawn.
The calculations are easiest to perform with a computer algebra system. In Mathematica a solution looks like this:
Compute the coefficients of the generating polynomials recursively:
prob[{j_, k___}, p_List, x_] := Block[{n = Total@p, q = p},
q[[j]] = 0;
prob[{k}, q, x] p[[j]] Subscript[x, j] / n] ;
prob[{}, _, _] := 1
Sum over all permutations to obtain the generating polynomials:
m[k_Integer, p_, x_] := m[k, p, x] =
Sum[prob[y, p, x], {y, Permutations[Range[Length[p]], {k}]}]
Specify the ideal generated by $x_1x_2\cdots x_n$ and the $x_i^2-x_i$:
ideal[p_, x_] :=
With[{n = Length[p]}, {Product[Subscript[x, i], {i, 1, n}],
Subscript[x, #]^2 - Subscript[x, #] & /@ Range[n]} // Flatten]
Multiply the generating polynomials, reducing modulo the ideal at each step:
Block[{weights = {4, 3, 1, 1, 1}, draws = {3, 3, 2}, x, i, gf},
i = ideal[weights, x];
gf = Fold[PolynomialMod[#1 m[#2, weights, x], i] &, 1, draws];
Set all $x_i$ to $1$:
gf /. (Subscript[x, #] -> 1 & /@ Range[Length[weights]])]
The answer takes $0.02$ seconds to compute and equals $\frac{74344589}{111132000}$. This means the chance of obtaining all five balls is $1 -\frac{74344589}{111132000} = \frac{36787411}{111132000} \approx 0.33102447.$
Discussion
Simulation
Whenever possible, it is a good idea to double-check complicated or tricky probability calculations with a simulation. Here is a quick simulation (for arbitrary weights and arbitrary numbers of draws) in R, as applied to the original problem:
weights <- c(4,3,1,1,1)
draws <- c(3,3,2)
n.iter <- 10^4
set.seed(17)
p <- mean(replicate(n.iter, {
length(table(unlist(sapply(draws,
function(d) sample.int(5, d, prob=weights))))) >= length(weights)
}))
cat("Estimate =", p,
"( Z =", (p - 36787411/111132000) / sqrt(p*(1-p)/n.iter), ")")
Estimate = 0.335 ( Z = 0.8422911 )
The output reports the proportion of iterations in which all the balls were observed, followed by a Z-statistic relative to the previously-computed answer. Because the Z-statistic is small, this is evidence that the computed answer is correct (up to a few times the standard error, which in this example is $0.005$).
Application to simpler problems
One case in which there is a formula (of sorts) is when there are no weights. The Mathematica code, applied to a weights vector of {1,1,1,1,1}, outputs $\frac{9}{20}$, yielding an answer of $1 - \frac{9}{20} = \frac{11}{20} = 0.55.$ Is it right?
The symmetry of this problem--all balls have equal probabilities--makes an independent calculation of the answer feasible. After drawing the first three balls, name them $x_1, x_2,$ and $x_3$, and let the remaining balls be $x_4$ and $x_5$. Replace all balls and draw another three. There are two relevant events to track: whether $3$ or $4$ balls have been observed. The chance of $3$ is the chance that the same $3$ balls were drawn a second time; this is $1 / \binom{5}{3} = 1/10.$ The chance of observing $4$ distinct balls is that two of the new balls were from $\{x_1, x_2, x_3\}$ and the third from $\{x_4, x_5\}$: this equals $\frac{\binom{3}{2}\binom{2}{1}}{\binom{5}{3}} = \frac{3\times 2}{10} = \frac{6}{10}.$
In each of these cases we compute the chance that the final draw of $2$ balls will not complete the set of five. If only $3$ balls have been observed, that chance is $9/10$. If $4$ balls have been observed, the chance is $\frac{\binom{4}{2}}{\binom{5}{2}} = \frac{6}{10}.$ Whence the chance of not completing the set is
$$\frac{1}{10} \frac{9}{10} + \frac{6}{10}\frac{6}{10} = \frac{9+36}{100} = \frac{45}{100} = \frac{9}{20},$$
exactly as reported by the Mathematica code.
Closer examination of the polynomials
How the code works is easy to see in this case, because the polynomials $m_k$ are multiples of the elementary symmetric functions:
$$10 m_2(x) = x_1 x_2+x_3 x_2+x_4 x_2+x_5 x_2+x_1 x_3+x_1 x_4+x_3 x_4+x_1 x_5+x_3 x_5+x_4 x_5$$
$$10 m_3(x) = x_1 x_2 x_3+x_1 x_4 x_3+x_2 x_4 x_3+\cdots+x_1 x_2 x_5+x_1 x_4 x_5+x_2 x_4 x_5.$$
These are constructed by taking a single product, such as $x_1x_2$, and adding to it all monomials attainable by permuting their subscripts (among the set $\{1,2,3,4,5\}$). In effect, each monomial has only one kind of term plus its permuted versions.
Let's compute their products modulo the ideal $\mathcal{I}$ generated by $\{x_1x_2x_3x_4x_5, x_1^2-x_1, x_2^2-x_2, x_3^2-x_3, x_4^2-x_4, x_5^2-x_5\}.$ To do so we only have to look at a limited number of products and then apply all permutations. For instance, in computing $(10 m_3(x))^2$ we find one term like $(x_1x_2x_3)^2 \equiv x_1x_2x_3$ mod $\mathcal{I}$--together with all its permutations--and terms that are congruent to $x_1x_2x_3x_4$ mod $\mathcal{I}$--together with all their permutations. It's straightforward to count that there are $12$ each of the latter, whence
$$(10 m_3(x))^2 \equiv\sum_{\sigma\in \mathcal{S}_5}\left(x_{\sigma(1)}x_{\sigma(2)}x_{\sigma(3)} + 12x_{\sigma(1)}x_{\sigma(2)}x_{\sigma(3)}x_{\sigma(4)}\right)\quad\text{mod }\mathcal{I}.$$
Multiplication by $10m_2(x)$ modulo $\mathcal{I}$ is similarly easy. The result is
$$(10 m_3(x))^2(10 m_2(x)) \equiv \sum_{\sigma\in \mathcal{S}_5}\left(3x_{\sigma(1)}x_{\sigma(2)}x_{\sigma(3)} + 84x_{\sigma(1)}x_{\sigma(2)}x_{\sigma(3)}x_{\sigma(4)}\right)\quad\text{mod }\mathcal{I}.$$
This says there are $\binom{5}{3}\times 3$ chances in $1000$ of observing exactly three of the balls after independent draws of $3$, $3$, and $2$ balls and $\binom{5}{4}\times 84$ chances in $1000$ of observing exactly four of the balls. Their sum is $10\times 3 + 5\times 84 = 450$ chances out of $1000$.
The general method used in the weighted-probability case performs a similar set of computations.
Alternative assumptions
For alternative drawing procedures (see the comment thread for the question), merely change the definition of prob, which describes the chances of drawing each sequence of balls. The rest of the algorithm remains unchanged.
