Motivated by applications to critical phenomena and open theoretical questions, we study conformal field theories with $O(m)times O(n)$ global symmetry in $d=3$ spacetime dimensions. We use both analytic and numerical bootstrap techniques. Using the analytic bootstrap, we calculate anomalous dimensions and OPE coefficients as power series in $varepsilon=4-d$ and in $1/n$, with a method that generalizes to arbitrary global symmetry. Whenever comparison is possible, our results agree with earlier results obtained with diagrammatic methods in the literature. Using the numerical bootstrap, we obtain a wide variety of operator dimension bounds, and we find several islands (isolated allowed regions) in parameter space for $O(2)times O(n)$ theories for various values of $n$. Some of these islands can be attributed to fixed points predicted by perturbative methods like the $varepsilon$ and large-$n$ expansions, while others appear to arise due to fixed points that have been claimed to exist in resummations of perturbative beta functions.