We used an optimal control method involving covariant control equations as optimality conditions, to command the actuators of robot manipulators. These form a coupled system of second order nonlinear ordinary differential equations when associated with the robot motion equations. By solving this system, the control action required to take the robot from an initial to a final state is optimized in a prescribed time. However, the target set of equations exhibited stiffness. Therefore, an adequate solution could only be found for short trajectory durations with readily available numerical methods. We examined a time discretization procedure based on cubic and quintic Hermite finite elements which exhibited superconvergence properties for interpolation. This motivated us to develop a time integration algorithm based on Hermite's technique, where motion and control equations were perturbed to solve the optimal control problem. The optimal motion of a robotic manipulator was simulated using this algorithm. Our method was compared with a commercial differential equations solver on the basis of specific indicators. It outperformed the commercial solver by effectively solving the stiff set of equations for longer trajectory durations, with the cubic elements performing better than the quintic ones in this sense. The convergence analysis of our method confirmed that the quintic elements are more precise at the cost of increased computational burden, but converge at a lower rate than expected. Controlled motion experiments on a robotic manipulator validated our methodology. Trajectories were smoothly tracked and results exposed further methodology improvements.