In this work, we aim at applying hybrid approaches for solving PDEs, combining machine learning and numerical methods. The goal is to build surrogate models that can accelerate computations and facilitate uncertainty analysis, while remaining compatible with stability and robustness constraints of numerical schemes. In this presentation, we propose an approach inspired by the iterative generalized polynomial chaos (i-gPC) method. This new approach involves progressively building and training a neural network with polynomial activations. The progressive construction eliminates the need to predefine the network architecture, thereby reducing hyperparameters and improving interpretability. This approach relies on optimization strategies based on the analysis of differential equations to ensure necessary theoretical properties. We demonstrate that this approach can tackle several limitations of i-gPC, such as instabilities and the lack of theoretical guarantees for strict error decay when a layer is added. The aim of this contribution is to present these approaches, study some numerical behaviors, and apply them to examples of physical interest.