Galactic systems are inherently multiphase. Understanding the roles of the various phases and their interactions is the next key step towards a more complete picture of galaxy evolution. A major challenge is that the transport and dynamics of cold clouds is governed by complex small scale processes. Large scale models thus require subgrid prescriptions in the form of models validated with small scale simulations. In this work, we explore using neural ordinary differential equations (NODEs), which embed a neural network, to more accurately model these clouds. We apply Symbolic Regression (SR) to potentially discover new insights into the physics of cloud-environment interactions. We test this on both generated mock data and actual simulation data. We also extend the model to include more than one neural term. We find that NODEs in tandem with SR can be used to enhance the accuracy and efficiency of subgrid models, and/or discover the underlying equations to improve generality and scientific understanding. We highlight the potential of this scientific machine learning approach as a natural extension to the traditional modelling paradigm, both for the development of semi-analytic models and for physically intepretable equation discovery in complex non-linear systems.