Langmuir probes are conceptually simple devices and have long been the primary diagnostic for plasmas, but interpreting the output is rarely straightforward, even in plasmas consisting of Maxwellian electrons and a single ion species. The problem is worse yet in more complicated plasmas where theoretical calculations are lacking. To help interpret probe data, this work presents a numerical scheme capable of modeling spherical and cylindrical probes in plasmas containing any number of charged species, positive or negative. The scheme is based on solving the derivative of Poisson's equation as an initial-value problem rather than solving Poisson's equation as a boundary-value problem. Because the new equation is linear in the electric field, stable solutions are possible. The same approach can be used with other nonlinear ordinary differential equations as well, but in most cases the initial conditions must be adjusted until the original boundary conditions are satisfied. Fortunately, this difficulty is easily circumvented for certain types of problems, including Langmuir probes. Accordingly, Langmuir probes can be modeled quickly and easily, even in plasmas containing multiple ion species and multi-temperature electrons. Examples are given for electropositive plasmas, electronegative plasmas, and two-temperature electrons.