Kohn-Sham (KS) inversion, in which the effective KS mean-field potential is found for a given density, provides insights into the nature of exact density functional theory (DFT) that can be exploited for the development of density functional approximations. Unfortunately, despite significant and sustained progress in both theory and software libraries, KS inversion remains rather difficult in practice, especially in finite basis sets. The present work presents a KS inversion method, dubbed the "Lieb-response"approach, that naturally works with existing Fock-matrix DFT infrastructure in finite basis sets, is numerically efficient, and directly provides meaningful matrix and energy quantities for pure-state and ensemble systems. Some additional work yields potential. It thus enables the routine inversion of even difficult KS systems, as illustrated in a variety of problems within this work, and provides outputs that can be used for embedding schemes or machine learning of density functional approximations. The effect of finite basis sets on KS inversion is also analyzed and investigated.