A simple, general, and numerically robust algorithm is presented for calculating the disjoining pressure and interaction free energy per unit area between two identically charged flat plates due to electrical double layer interactions according to the nonlinear Poisson-Boltzmann theory. The result is applicable to electrolytes with any number of ionic species having any combination of valencies as well as to constant potential, constant charge, or charge regulation boundary conditions on the plates. The algorithm is very simple to implement on commonly available numerical software environments and is therefore particularly suitable for use in data analysis.