When performing long-time predictions or other extrapolation tasks, knowing an accurate ODE or PDE model can be crucial, especially when there is a limited amount of data. However, when the data is scarce, knowing which model is best suited to our data is not an easy task. In this work, we propose a method to asses the bias of a given model, using the available noisy data, and to learn a correction term that is then added to the model, making it better adapted to our measurements. It is a nonlinear operator that maps the model solution and the model parameters to a correction function. For example, this correction term could act on poorly chosen boundary conditions or can be a volumetric term. Once this term is learned, the corrected model can be solved directly for extrapolation tasks. If the model is biased, the model parameters associated with the biased model are typically unsuitable for the corrected model. Therefore, when learning the correction operator we must simultaneously learn information on the model parameters. In this work, we choose a Bayesian estimation method, learning a posterior distribution on the model parameters given the data and the corrective term. We use the Stein Variational Gradient Descent (SVGD) to sample from the posterior distribution. Due to high complexity when a large number of samples are required, we propose an optimally relaxed variant of the SVGD, which accelerates the sampling by up to 2 orders of magnitude. Our method is agnostic of the form given to the correction term: it can be expressed using basis functions, tensors, or neural networks, with special interest in neural operator architectures.