The paper presents a gradient-based numerical algorithm for optimal control of nonlinear multivariable systems with control and state vectors constraints. The algorithm has a backward-in-time recurrent structure similar to the backpropagation-through-time (BPTT) algorithm, which is mostly used as a learning algorithm for dynamic neural networks. This paper presents an enhancement of the basic optimization algorithm. Our enhanced algorithm uses high-order Adams time-discretization schemes instead of the basic Euler discretization method, and a numerical calculation of Jacobians as an alternative to analytical Jacobians. Two examples are considered to illustrate the algorithm and its performance. The first example is that of a tubular reactor, for which an analytical solution is available, which can be readily used for validation of our approach. The second example is related to controlling vehicle dynamics based on a realistic high order model.