A variable-grid, high-order finite difference (FD) method is applied to the modeling of mantle convection in both two- and three-dimensional geometries. The algorithm combines extreme simplicity in programming with a very high degree of accuracy. Memory requirements are low and grow almost linearly with the total number of grid points in three dimensions, regardless of the increase in grid points in the vertical direction. Higher-order methods, such as eighth order, yield significantly better results than a second-order method for the same grid size, with only a modest increase in memory requirements. This is particularly important for high Rayleigh number convection, where the large number of grid points required to obtain an accurate enough solution with second-order schemes would make the computation extremely costly. The small-scale features in the hard-turbulent regime under high-Rayleigh number situations can greatly stress low-order methods, and in these situations a high-order method is definitely needed. We have numerically simulated three-dimensional time-dependent convection for constant properties up toRa= 108, using an eighth-order FD scheme. Both purely base-heated and partially internally heated situations have been considered. The hot plumes are broader near the surface with internal heating. Detailed studies of the three-dimensional constant viscosity plumes indicate that no small-scale circulation takes place in the ascending plume heads regardless of the heating configuration in accordance with predictions from boundary layer theory.