One of the efficient and powerful schemes to solve linear and nonlinear equations is homotopy analysis method (HAM). In this work, we obtain the approximate solution of a system of partial differential equations (PDEs) by means of HAM. For this purpose, we develop the concept of HAM for a system of PDEs as a matrix form. Then, we prove the convergence theorem and apply the proposed method to find the approximate solution of some systems of PDEs. Also, we show the region of convergence by plotting the H-surface.