The self- and mutual-avoiding effects of two star polymers in a good solvent are studied by means of a simple Monte Carlo sampling technique on a cubic lattice, using an enrichment algorithm. The total number of configurations is determined as a function of the distance between the two star polymers. The radius of gyration of a single star polymer, the second virial coefficient, the effective interstar potential, and the penetration function are evaluated for 3- to 6-arm star polymers. Considerable deviation from the result of the previous naive first-order epsilon-expansion is observed in the penetration function. We also discuss the value of the penetration function in the infinite-arm limit, according to the cone picture.