This research is dedicated to the discrete-time tracking control of nonlinear multi-input multi-output (MIMO) Newtonian mechanics systems with an equal number of inputs and outputs, in the presence of input saturation and high-order dynamics of the actuation system. Firstly, the continuous-time model of the plant and actuation system has been transformed into the discrete-time domain using Adams–Bashforth and Adams–Moulton discretization techniques. In the next step, considering saturation boundaries, the modified position and velocity references are derived using the discrete-time reaching law-based control approach. These modified references are obtained through an optimization algorithm that calculates the closest achievable position and velocity trajectories with respect to the desired reference trajectory. Modifying desired trajectories in MIMO systems results in the formation of a system of linear matrix inequalities (LMIs). Solving such systems of inequalities requires the use of complex methods, and this difficulty increases with the increase in the system’s degrees of freedom (DOFs). The proposed strategy offers a more straightforward approach to solving these inequalities. Once the modified references are calculated, the corresponding position and velocity control signals can be obtained using the mentioned reaching law strategies. The final step of the proposed algorithm is to combine the resulting position and velocity control signals using a Gaussian-weighted averaging method to form the main control signal. A desirable performance for the control system can be achieved by appropriately adjusting the standard deviations of the Gaussian weighing functions. To evaluate the performance characteristics of the proposed method, it has been applied in the tracking control procedure of a nonlinear 3-DOF mechanical system. The numerical simulation results demonstrate the effectiveness of the method.