The bulk kinematics and thermodynamics of hot supernovae-driven galactic winds is critically dependent on both the amount of swept up cool clouds and non-spherical collimated flow geometry. However, accurately parameterizing these physics is difficult because their functional forms are often unknown, and because the coupled non-linear flow equations contain singularities. We show that deep neural networks embedded as individual terms in the governing coupled ordinary differential equations (ODEs) can robustly discover both of these physics, without any prior knowledge of the true function structure, as a supervised learning task. We optimize a loss function based on the Mach number, rather than the explicitly solved-for 3 conserved variables, and apply a penalty term towards near-diverging solutions. The same neural network architecture is used for learning both the hidden mass-loading and surface area expansion rates. This work further highlights the feasibility of neural ODEs as a promising discovery tool with mechanistic interpretability for non-linear inverse problems.