Many scientific and engineering challenges require tools capable of discovering unknown, intricate patterns and relationships in the form of symbolic models. These symbolic models are favored for their interpretability and explanatory power that are vital for understanding complex systems and making well-informed decisions. Although modern large deep learning models, particularly transformer-based architectures, demonstrate impressive performance in symbolic regression tasks, they heavily depend on the training data and are therefore poorly transferable to new domains. In contrast, traditional online optimization approaches, such as genetic programming algorithms, are independent of training data and enable the search for out-of-distribution models. However, convergence to the optimal solution with genetic programming algorithms is not guaranteed. Meanwhile, methods based on linear regression, though computationally demanding, are constrained by the proper choice of basis expressions, which may require domain specific knowledge. In this work, we propose a novel algorithm for symbolic regression and model discovery. Leveraging the concept of Neural Operators (NOs), known for their capability to approximate non-linear operators, we first train Fourier Neural Operators (FNOs) to approximate a set of library functions. Subsequently, we construct a neural network from these trained operators. To identify the optimal symbolic expression, we search for a sparse combination of weights that best fits the data through a backward propagation. The weights obtained from optimization, along with the network structure, can then be interpreted as a symbolic expression. We evaluate the proposed approach on well-known symbolic regression benchmarks and compare its performance with the current state-of-the-art methods.