Discovering the true underlying pharmacological processes that evolve over time is of significant interest to pharmacometricians and healthcare professionals. Accurate and interpretable modeling of these processes is crucial for understanding disease-drug interactions, optimizing medication dosing, and minimizing adverse health effects, such as in chemotherapy where precise dosing can maximize cancer cell eradication while minimizing patient harm. Current models, often developed by human experts, are limited by their high cost, lack of scalability, and restriction to existing human knowledge. In this paper, we present the Data-Driven Discovery (D3) framework, a novel approach leveraging Large Language Models (LLMs) to iteratively discover and refine interpretable models of pharmacological dynamics. Unlike traditional methods, D3 enables the LLM to propose, acquire and integrate new features, validate, and compare pharmacological dynamical systems models, thereby uncovering new insights into pharmacokinetic processes. We demonstrate that D3 can identify well-fitting, interpretable models across diverse pharmacokinetic datasets. Specifically, experiments on a real pharmacokinetic Warfarin dataset reveal that D3 uncovers a new plausible pharmacokinetic model that outperforms existing literature, highlighting its potential for precision dosing in clinical applications.