Discovering mathematical models that characterize the observed behavior of dynamical systems remains a major challenge, especially for systems in a chaotic regime, due to their sensitive dependence on initial conditions and the complex non-linear interactions within the system. The challenge is even greater when the physics underlying such systems is not yet understood, and scientific inquiry must solely rely on empirical data. Despite advancements such as the Sparse Identification of Nonlinear Dynamics (SINDy) algorithm, which has shown success in identifying chaotic systems with reasonable accuracy, challenges remain in dealing with noise, sparse data, and the need for models that generalize well across different chaotic systems. Driven by the need to fill this gap, we develop a framework named AI-Lorenz that learns mathematical expressions modeling complex dynamical behaviors by identifying differential equations from noisy and sparse observable data. We train a physics-informed neural network (PINN) with the eXtreme Theory of Functional Connections (X-TFC) algorithm by using data and known physics (when available) to learn the dynamics of a system, its rate of change in time, and missing model terms, which are used as input for a symbolic regression algorithm (PySR) to autonomously distill the explicit mathematical terms. This, in turn, enables us to predict the future evolution of the dynamical behavior. The performance of this framework is validated by recovering the right-hand sides and unknown terms of certain complex, chaotic systems, such as the well-known Lorenz system, a six-dimensional hyperchaotic system, the non-autonomous Sprott chaotic system, and the slow-fast Duffing system, and comparing them with their known analytical expressions and state-of-the-art regression and system identification methods.