We combined multiple final states ensuing from non-resonant Higgs pair production by optimising over each final state using multivariate techniques and put bounds on the self-coupling of Higgs by employing the log-likelihood confidence level hypothesis test. Further, we calculate the production cross-section limits from multifarious heavy Higgs decay channels, including resonant di-Higgs production, in a model-independent way and draw its impact on the MSSM parameter space at the HL-LHC.