We present a robust method for computing the reflection of arbitrarily shaped and sized beams from finite thickness photonic crystals. The method is based on dividing the incident beam into plane waves, each of which can be solved individually using Bloch periodic boundary conditions. This procedure allows us to take a full advantage of the crystal symmetry and also leads to a linear scaling of the computation time with respect to the number of plane waves needed to expand the incident beam. The algorithm for computing the reflection of an individual plane wave is also reviewed. Finally, we find an excellent agreement between the computational results and measurement data obtained from opals that are synthesized using polystyrene and poly(methyl methacrylate) microspheres.