The role of higher order coupling of surface vibrations to the relative motion in heavy-ion fusion reactions at near-barrier energies is investigated. The coupled channels equations are solved to all orders and also in the linear and the quadratic coupling approximations. Taking ${}^{64}$Ni + ${}^{92,96}$Zr reactions as examples, it is shown that all order couplings lead to considerably improved agreement with the experimentally measured fusion cross sections and average angular momenta of the compound nucleus for such heavy, nearly symmetric systems. The importance of higher order coupling is also examined for asymmetric systems like ${}^{16}$O + ${}^{112}$Cd, ${}^{144}$Sm, for which previous calculations of the fusion cross section seemed to indicate that the linear coupling approximation was adequate. It is shown that the shape of the barrier distributions and the energy dependence of the average angular momentum can change significantly when the higher order couplings are included, even for systems where measured fusion cross sections may seem to be well reproduced by the linear coupling approximation.