Abstract Constitutive modeling is the cornerstone of computational and structural mechanics. In a finite element analysis, the constitutive model is encoded in the material subroutine, a function that maps local strains onto stresses. This function is called within every finite element, at each integration point, within every time step, at each Newton iteration. Today’s finite element packages offer large libraries of material subroutines to choose from. However, the scientific criteria for model selection remain highly subjective and prone to user bias. Here we fully automate the process of model selection, autonomously discover the best model and parameters from experimental data, encode all possible discoverable models into a single material subroutine, and seamlessly integrate this universal material subroutine into a finite element analysis. We prototype this strategy for tension, compression, and shear data from human brain tissue and perform a hyperelastic model discovery from twelve possible terms. These terms feature the first and second invariants, raised to the first and second powers, embedded in the identity, exponential, and logarithmic functions, generating 2 2×2×3 = 4096 models in total. We demonstrate how to integrate these models into a single universal material subroutine that features the classical neo Hooke, Blatz Ko, Mooney Rivlin, Demiray, Gent, and Holzapfel models as special cases. Finite element simulations with our universal material subroutine show that it specializes well to these widely used models, generalizes well to newly discovered models, and agrees excellently with both experimental data and previous simulations. It also performs well within realistic finite element simulations and accurately predicts stress concentrations in the human brain for six different head impact scenarios. We anticipate that integrating automated model discovery into a universal material subroutine will generalize naturally to more complex anisotropic, compressible, and inelastic materials and to other nonlinear finite element platforms. Replacing dozens of individual material subroutines by a single universal material subroutine that is populated directly via automated model discovery—entirely without human interaction—makes finite element analyses more accessible, more robust, and less vulnerable to human error. This could forever change how we simulate materials and structures.