Symbolic regression is a machine learning technique that can learn the governing formulas from data and thus has the potential to transform scientific discovery. However, symbolic regression is still limited in the complexity of the systems that it can analyze. Deep learning on the other hand has transformed machine learning in its ability to analyze extremely complex and high-dimensional datasets. Here we develop a method that uses neural networks to extend symbolic regression to parametric systems where some coefficient may vary as a function of time but the underlying governing equation remains constant. We demonstrate our method on various analytic expressions and PDEs with varying coefficients and show that it extrapolate well outside of the training domain. The neural network-based architecture can also integrate with other deep learning architectures so that it can analyze high-dimensional data while being trained end-to-end in a single step. To this end we integrate our architecture with convolutional neural networks and train the system end-to-end to discover various physical quantities from 1D images of spring systems where the spring constant may vary.