We consider the Stefan problem in a domain Ω ⊆ R^d, (d=1,2,3) ∂tu(·; μ)-∆β(u(·; μ); μ) = f (·; μ) in Ω × (0, T ) supplemented with initial and boundary conditions. Here β(·; μ) : R → R is a (non strictly) increasing Lipschitz function. We assume that the initial and boundary conditions, as well as f and β depend on some parameter μ ∈ P ⊆ R^p, p > 1. Such problem is typically used to describe solidification, arising e.g. during laser additive processes. In this context, the function β is given by the properties of the material, whereas the laser heat source is modelled by the source term f . The Stefan problem being nonlinear, classical numerical methods used to approximate its solution are time consuming. This becomes an issue in the many query context, e.g. if an inverse problem has to be solved, or if the solution corresponding to a new parameter μ is needed in real time. To address this problem, we aim to build a neural network that approximates the map (x, t; μ) → u(x, t; μ). For this purpose, we consider a data-driven neural network trained with data coming from finite element simulations u_h(·; μi), i = 1, . . . N , where h denotes the numerical parameters (mesh size and time step). We denote the neural network approximation by u_N . To ensure that the network gives a reliable approximation of u, we want to ensure that the error between u and u_N is below some preset tolerance. For this purpose, we decompose the overall error in a given norm || · || over Ω × (0, T ) × P as ||u − u_N || ≤ ||u − u_h|| + ||u_h − u_N ||. The first term corresponds to the finite element error, that is embedded in the training data. To control it, we use a mesh adaptation algorithm relying on an a posteriori error indicator. We propose an adaptive algorithm that increases the size of the training set to control the error of the neural network. Once a (sufficiently accurate) neural network has been obtained, it can be used to efficiently compute the solution corresponding to a new parameter μ, or to solve an inverse problem. Numerical results will be presented for a 2D model problem.