We propose a new algorithm for computing the luminosity distance in the flat universe with a cosmological constant based on Shchigolevs homotopy perturbation method, where the optimization idea is applied to prevent the arbitrariness of initial value choice in Shchigolevs homotopy. Compared with the some existing numerical methods, the result of numerical simulation shows that our algorithm is a very promising and powerful technique for computing the luminosity distance, which has obvious advantages in computational accuracy,computing efficiency and robustness for a given {Omega_m}.