This is a Monte Carlo program to estimate the probabilities in problem 18 of the Fall 2011 MAT 201 course.

This chooses three distinct integers in the range 1, 2, 3, 4, 5.  I will use this to pick three balls at random.  I imagine putting them in a long line, with the r reds first and the w whites following.  The three integers are ordered with 1<=  i < j < k <= n = r+w.
if r < i, in which case r < i, j, k and all the balls are white, then throw away the white balls.
else if r < j, in which case i <= r < j < k and one is red and 2 are white, then throw away the red ball.
else if r < k, in which case i < j <= r < k and two are red and one is white, then throw away the white ball.
else if k <= r, in which cas i < j < k <= r and all the balls are red, then throw awway the 3 red balls.

RandomSample[Range[5],3]

pull3balls[rwList_] := Module[{r, w, n,newR, newW, i, j, k},
{r, w} = rwList;(* number of red and white balls*)
n = r+w;
If[n <3, Return[rwList]];
{i,j,k} = Sort@RandomSample[Range[n], 3];
(*Print[{i,j,k}];*)
Which[r < i, (*all balls are white*) newR = r; newW= w-3,
r < j (* one red and 2 white *), newR = r-1; newW = w,
r < k (* two red and 1 white *), newR = r; newW = w-1,
k ≤ r (* all the balls are red *), newR = r-3; newW = w];
(*output new number of balls *)
{newR, newW} (* this is the output of the function *)
]

I put in the Print[{i, j, k}]; command and ran this several times to check that it works.

pull3balls[{2,0}]

playOneGame[rwListInput_] := Module[{keepGoing, rwList,newRwList},
rwList =rwListInput;
keepGoing = True;
While[keepGoing,
(*Print[rwList];*)
newRwList =  pull3balls[rwList] ;
keepGoing = (newRwList ≠ rwList);
rwList = newRwList];
rwList (*output*)
]

I did several of these with the print statement not commented out.

playOneGame[{5,5}]

This next cell is the  actual run.

nTrials = 100;
rwList = {2010,2010};
Print["nTrials = ", nTrials, ", and initial rwList = ", rwList];
n00 = 0;
n11 = 0;
n01 = 0;
n10 = 0;
n20 = 0;
n02 = 0;
{sec, dummy} =Timing[
Do[result =playOneGame[rwList];
Which[
result == {0,0}, n00 ++,
result == {1,1}, n11 ++,
result == {0,1}, n01 ++,
result == {1,0}, n10 ++,
result == {2,0}, n20 ++,
result == {0,2}, n02 ++];
(*Print[result]*), {nTrials}]];
nt = N[nTrials];
Print["P00 = ", n00/nt, ", P11 = ", n11/nt, ", P01 = ", n01/nt, ", P10 = ", n10/nt, ", P20 = ", n20/nt, ", P02 = ", n02/nt];
Print["This run took ", sec, " seconds."];

Here is a rather short run, with 2010 balls of each color.  It doesn’t seem different from the 2011 ball problem.

This is the result of an overnight run.  It took about 5 and 1/4 hours.

19016/3600.

A very rough estimtmate of the error is the difference between P01 and P10, and between P20 and P02.  It is also about 1/sqrt [nTrials]. They all give an error in the third decimal place.

n10/nt - n01/nt
n20/nt - n02/nt

1/Sqrt[100000.]

122/1000

825/10000

The estimates of the probabilities shown retain the full resolution of the results:

Print[{n00, n11, n10, n10, n20, n02}]