  A completely integrable nonlinear partial differential equation (PDE) can be
associated with a system of linear PDEs in an auxiliary function whose
compatibility requires that the original PDE is satisfied. This associated
system is called a Lax pair. Two equivalent representations are presented. The
first uses a pair of differential operators which leads to a higher order
linear system for the auxiliary function. The second uses a pair of matrices
which leads to a first-order linear system. In this paper we present a method,
which is easily implemented in Maple or Mathematica, to compute an operator Lax
pair for a set of PDEs.
  In the operator representation, the determining equations for the Lax pair
split into a set of kinematic constraints which are independent of the original
equation and a set of dynamical equations which do depend on it. The kinematic
constraints can be solved generically. We assume that the operators have a
scaling symmetry. The dynamical equations are then reduced to a set of
nonlinear algebraic equations.
  This approach is illustrated with well-known examples from soliton theory. In
particular, it is applied to a three parameter class of fifth-order KdV-like
evolution equations which includes the Lax fifth-order KdV, Sawada-Kotera and
Kaup-Kuperschmidt equations. A second Lax pair was found for the Sawada--Kotera
equation.
