Despite the rich amount of information that can be extracted by the 2DES signals, there is only a limited amount of theoretical work on simulating 2DES for exciton polariton. In our previous work , we developed a theoretical approach for simulating linear and nonlinear spectroscopy of exciton-polaritons within the semiclassical partially linearized path-integral framework, but it was limited to N = 1 molecule. This is because, with an increasing number of molecules N , the number of states within the second excitation manifold expands quadratically in order of O(N 2 ) and the Hamiltonian size scales as O(N 4 ). Due to this quartic scaling, the direct computation of dynamical properties for an N -molecule polaritonic system becomes very expensive even for a system as small as N = 10. ical development in Ref. 8 and in paper-I 1 to simulate 2DES of N -molecules polariton systems, while introducing new algorithms for the focusing step in L-PLDM 2DES calculations. These new developments reduce the cost of the calculations to the order O(N 2 ) without the necessity for generating any O(N 4 ) matrix. We further show that our algorithm can be efficiently vectorized using GPU to reduce the scaling of 2DES propagation time R from O(T 3 ) to O(T 0 ), where T is the number of propagation steps for each laser in the simulation. This allows us to efficiently perform accurate 2DES calculations, with the help of GPU, for polariton systems with N ≈ O(10 2 ) where the total states can be of order O ). We demonstrate the 2DES for up to N = 25 molecules coupled to the cavity at t 2 =0 fs and also demonstrate the signal breakdown of 2DES for N = 5 molecule polariton system at different population time t 2 . Finally, we analyzed the features of different signals from both rephasing and non-rephasing contributions while also further decomposing them as ground-state bleaching (GSB), excited state emission (ESA), and stimulated emission (SE) pathways.