Invariant solutions of the Navier-Stokes equations play an important role in the spatiotemporally chaotic dynamics of turbulent shear flows. Despite the significance of these solutions, their identification remains a computational challenge, rendering many solutions inaccessible and thus hindering progress towards a dynamical description of turbulence in terms of invariant solutions. We compute equilibria of three-dimensional wall-bounded shear flows using an adjoint-based matrix-free variational approach. To address the challenge of computing pressure in the presence of solid walls, we develop a formulation that circumvents the explicit construction of pressure and instead employs the influence matrix method. Together with a data-driven convergence acceleration technique based on dynamic mode decomposition, this yields a practically feasible alternative to state-of-the-art Newton methods for converging equilibrium solutions. We compute multiple equilibria of plane Couette flow starting from inaccurate guesses extracted from a turbulent time series. The variational method outperforms Newton(-hookstep) iterations in converging successfully from poor initial guesses, suggesting a larger convergence radius.