  A recursion operator is an integro-differential operator which maps a
generalized symmetry of a nonlinear PDE to a new symmetry. Therefore, the
existence of a recursion operator guarantees that the PDE has infinitely many
higher-order symmetries, which is a key feature of complete integrability.
Completely integrable nonlinear PDEs have a bi-Hamiltonian structure and a Lax
pair; they can also be solved with the inverse scattering transform and admit
soliton solutions of any order.
  A straightforward method for the symbolic computation of polynomial recursion
operators of nonlinear PDEs in (1+1) dimensions is presented. Based on
conserved densities and generalized symmetries, a candidate recursion operator
is built from a linear combination of scaling invariant terms with undetermined
coefficients. The candidate recursion operator is substituted into its defining
equation and the resulting linear system for the undetermined coefficients is
solved.
  The method is algorithmic and is implemented in Mathematica. The resulting
symbolic package PDERecursionOperator.m can be used to test the complete
integrability of polynomial PDEs that can be written as nonlinear evolution
equations. With PDERecursionOperator.m, recursion operators were obtained for
several well-known nonlinear PDEs from mathematical physics and soliton theory.
