Partial differential equations are transformed into ordinary differential equations, and a variational formulation is then established. The trial function is chosen using Jacobi-elliptic function with some unknown parameters similar to the exp-function method. Various approximate solitary solutions are obtained when making the obtained variational principle stationary with respect to each unknown parameter in the trial-function. The coupled Zakharov-Kuznetsov equations are used as an example to elucidate the solution procedure.