We apply the notion of discrete supersymmetry based on matrix factorization to quantum systems consisting of coupled bosonic oscillators to construct isospectral bosonic quantum networks. By using the algebra that arises due to the indistinguishability of bosonic particles, we write down the Schr dinger equations for these oscillators in the different boson-number sectors. By doing so, we obtain, for every partner quantum network, a system of coupled differential equations that can be emulated by classical light propagation in optical waveguide arrays. This mathematical scheme allows us to build quasi-two-dimensional optical arrays that are either isospectral or share only a subset of their spectrum after deliberately omitting some chosen eigenstates from the spectrum. As an example, we use this technique (which we call bosonic discrete supersymmetry or BD-SUSY) to design two optical, silica-based waveguide arrays consisting of six and three elements, respectively, with overlapping eigenspectrum.