Discrete variational methods show excellent performance in numerical simulations of different mechanical systems. In this paper, we introduce an iterative procedure for the solution of discrete variational equations for boundary value problems. More concretely, we explore a parallelization strategy that leverages the capabilities of multicore CPUs and GPUs (graphics cards). We study this parallel method for higher-order Lagrangian systems, which appear in fully-actuated problems and beyond. The most important part of the paper is devoted to a precise study of different convergence conditions for these methods. We illustrate their excellent behavior in some interesting examples, namely Zermelo's navigation problem, a fuel-optimal navigation problem, interpolation problems or in a fuel optimization problem for a controlled 4-body problem in astrodynamics showing the potential of our method.