We represent low dimensional quantum mechanical Hamiltonians by moderately sized finite matrices that reproduce the lowest O(10) bound-state energies and wave functions to machine precision. The method extends also to Hamiltonians that are neither Hermitian nor PT symmetric and thus allows one to investigate whether or not the spectra in such cases are still real. Furthermore, the approach is especially useful for problems in which a position-dependent mass is adopted, for example in effective-mass models in solid-state physics or in the approximate treatment of coupled nuclear motion in molecular physics or quantum chemistry. The performance of the algorithm is demonstrated by considering the inversion motion of different isotopes of ammonia molecules within a position-dependent mass model and some other examples of one- and two-dimensional Hamiltonians that allow for the comparison to analytical or numerical results in the literature.