The two-dimensional (2-D) parabolic equation (PE) is widely used for making radiowave propagation predictions in the troposphere. The effects of transverse terrain gradients, propagation around the sides of obstacles, and scattering from large obstacles to the side of the great circle path are not modeled, leading to prediction errors in many situations. In this paper, these errors are addressed by extending the 2-D PE to three dimensions. This changes the matrix form of the PE making it difficult to solve. A novel iterative solver technique, which is highly efficient and guaranteed to converge, is being presented. In order to confine the domain of computation, a three-dimensional (3-D) rectangular box is placed around the region of interest. A new second-order nonreflecting boundary condition is imposed on the surface of this box and its angular validity is established. The boundary condition is shown to keep unwanted fictitious reflections to an acceptable level in the domain of interest. The terrain boundary conditions for this 3-D PE method are developed and an original technique for incorporating them into the matrix form of the iterative solver is described. This is done using the concept of virtual field points below the ground. The prediction accuracy of the 3-D PE in comparison to the 2-D PE is tested both against indoor scaled frequency measurements and very high frequency (VHP) field trials.