We present a computational framework for canonical quantization in arbitrary inhomogeneous dielectric media by incorporating quantum electromagnetic effects into complex solutions of quantum Maxwell's equations. To do so, the proposed algorithm integrates and performs (1) numerical computation of normal modes and (2) evaluation of arbitrary products of ladder operators acting on multimode Fock states. The former is associated with Hermitian-Helmholtz linear systems using finite-element or finite-difference methods; consequently, the complete set of numerical normal modes diagonalizes the Hamiltonian operators up to floating-point precision. Its Hermiticity is retained, allowing its quantization. Then, we perform quantum numerical simulations of two-photon interference occurring in a 50:50 beam splitter to observe the Hong-Ou-Mandel effect. Our prototype model is useful for numerical analyses on various narrow-band quantum-optical multiphoton systems such as quantum metasurfaces, quantum-optical filters, and quantum electrodynamics in open optical cavities.