In this article we have introduced the mathematical theory of multistep optimal homotopy asymptotic method (MOHAM). The proposed method is implemented to different models having system of partial differential equations (PDEs). The results obtained by proposed method are compared with Homotopy Analysis Method (HAM) and closed form solutions. The comparisons of these results show that MOHAM is simpler in applicability, effective, explicit, control the convergence through optimal constants, involve less computational work. The MOHAM is independent of the assumption of initial conditions and small parameters like Homotopy Perturbation Method (HPM), HAM, Variational Iteration Method (VIM), Adomian Decomposition Method (ADM) and Perturbation Method (PM).