Conventional methods for numerical upper bound limit analysis utilize a predefined mesh of elements with nodes at fixed coordinates, and there the objective is to find the kinematically admissible velocity field associated with the mesh that minimizes the limit load. In this paper, a numerical method for optimizing the applied load with respect to both the velocity field and nodal coordinates is presented. The formulation is for the classical case of rigid elements (blocks) separated by velocity discontinuities, considering plane strain and weightless material obeying the Mohr-Coulomb yield condition. The central idea of the approach is to successively perturb velocities and nodal coordinates starting from an initial solution corresponding to a predefined mesh, where the equations governing the perturbed velocities and nodal coordinates are derived through linearization of exact expressions of the limit load and the conditions imposed on velocity jumps by the flow rule. Despite approximations made in perturbation steps, the method furnishes a final solution which satisfies the flow rule and energy balance exactly (within numerical tolerances) and therefore provides a rigorous upper bound on the true collapse load. The effectiveness of the method is demonstrated by comparing numerical results with the analytical solution for a benchmark problem in soil mechanics.