The finite difference method (FDM) is used to compute electrostatic potential distributions in photomultipliers. New Fortran packages, using both successive over-relaxation (SOR) and incomplete Choleski conjugate gradient (ICCG) techniques, have been developed for solving the finite difference equations. The effects of the mesh size on the accuracy of the results and the difference between the two methods are highlighted. The electron trajectories are computed by direct ray tracing with a power series method. The software can handle electron transparent grids as well as dynodes. It has been used to characterize the performance of many photomultiplier tubes. The results agree very well with experimental measurements. In addition, the advantages of using three-dimensional field computation in the design of certain PM tubes are illustrated.