The paper addresses the problem of accelerating predictive control of non-linear system models augmented with Gaussian processes (GP-MPC). Due to the non-linear and stochastic prediction model, predictive control of GP-based models requires to solve a stochastic optimization problem. Different model simplification methods have to be applied to reformulate this problem to a deterministic, non-linear optimization task that can be handled by a numerical solver. As these problems are still complex, especially with exact moment calculations, real-time implementation of GP-MPC is extremely challenging. The existing solutions accelerate the computations at the solver level by linearizing the non-linear optimization problem and applying sequential convexification. In contrast, this paper proposes a novel GP-MPC solution approach that without linearization formulates a series of surrogate quadratic programs (QP-s) to iteratively obtain the solution of the original non-linear optimization problem. The first step is embedding the non-linear mean-variance dynamics of the GP-MPC prediction model in a linear parameter-varying (LPV) structure and rewriting the constraints in parameter-varying form. By fixing the scheduling trajectory at a known variation (based on previously computed or initial state-input trajectories), optimization of the input sequence for the remaining varying linear model reduces to a linearly constrained quadratic program. After solving the QP, the non-linear prediction model is simulated for the new control input sequence and new scheduling trajectories are updated. The procedure is iterated until the convergence of the scheduling, that is, the solution of the QP converges to the solution of the original non-linear optimization problem. By designing a reference tracking controller for a 4DOF robot arm, we illustrate that the convergence is remarkably fast and the approach is computationally advantageous compared to current solutions. The proposed method enables the application of GP-MPC algorithms even with exact moment matching on fast dynamical systems and requires only a QP solver.