Summary: This work introduces a surrogate modeling strategy, based on diffusion maps manifold learning and artificial neural networks. On this basis, a numerical procedure is developed for cost-efficient predictions of a complex system’s response modeled by parametrized partial differential equations. The idea is to utilize a collection of solution snapshots, obtained by solving the partial differential equation for a small number of parameter values, in order to establish an efficient yet accurate mapping from the problem’s parametric space to its solution space. In common practice, solving the partial differential equations in the framework of the finite element method leads to high-dimensional data sets, which are a major challenge for machine learning algorithms to handle (curse of dimensionality). To overcome this problem, the proposed method exploits the dimensionality reduction properties of the diffusion maps algorithm in order to obtain a meaningful low-dimensional representation of the solution data set. With this approach, a reduced set of ‘hyperparameters’ is obtained, namely, the diffusion map coordinates, that characterize the high-dimensional solution vectors. Using this reduced representation, a feed-forward neural network is efficiently trained that maps the problem’s parameter values to their images in the low-dimensional diffusion maps space. This approach offers two advantages. The obvious one is that it reduces the computational resources and CPU time needed to train the neural network, which can be prohibitive for high dimensional problems. The second advantage is that training the neural network on the diffusion map coordinates, essentially translates to using the diffusion maps distance in the MSE loss function of the network. Compared to the Euclidean distance in the ambient space, the diffusion distance gives a better approximation of the distance between two points on the solution manifold, which leads to a more accurate network. Lastly, a mapping is developed based on the Laplacian pyramids scheme in order to convert points from the diffusion maps space back to the solution space. The composition of the neural network with the Laplacian pyramid scheme, substitutes the process of solving the partial differential equation under consideration, and is capable of providing the additional system solutions at a very low cost and high accuracy.