In this paper, the Optimal Homotopy Asymptotic Method is extended to derive the approximate solutions of fractional order two-dimensional partial differential equations. The fractional order Zakharov–Kuznetsov equation is solved as a test example, while the time fractional derivatives are described in the Caputo sense. The solutions of the problem are computed in the form of rapidly convergent series with easily calculable components using Mathematica. Reliability of the proposed method is given by comparison with other methods in the literature. The obtained results showed that the method is powerful and efficient for determination of solution of higher-dimensional fractional order partial differential equations.