Multi-channel optical systems such as disordered media, metasurfaces, multi-mode fibers, and photonic circuits are of both fundamental and technological importance. However, their optical responses are difficult to describe in the presence of multiple scattering, strong coupling, and multipath interference. Accurate modeling requires many large-scale simulations associated with the many input states, which currently take enormous computing resources. We propose an approach where all simulations are solved jointly and efficiently by augmenting the Maxwell operator with the source and the projection profiles, followed by a single partial factorization. This method applies to any linear system with any type of inputs, and benchmarks show it uses less memory and is 1,000 to 30,000,000 times faster than existing methods. It opens the door to previously inaccessible studies across disciplines involving multi-channel wave transport. As pilot examples, we realize, for the first time, full-wave simulations of entangled-photon backscattering from disorder and all-angle characterizations of aperiodic metasurfaces that are thousands of wavelengths wide.