It is now possible to compute linear in mass-ratio terms in the post-Newtonian (PN) expansion for compact binaries to very high orders using linear black hole perturbation theory applied to various invariants. For instance, a computation of the redshift invariant of a point particle in a circular orbit about a black hole in linear perturbation theory gives the linear-in-mass-ratio portion of the binding energy of a circular binary with arbitrary mass ratio. This binding energy, in turn, encodes the system's conservative dynamics. We give a method for extracting the analytic forms of these post-Newtonian coefficients from high-accuracy numerical data using experimental mathematics techniques, notably an integer relation algorithm. Such methods should be particularly important when the calculations progress to the considerably more difficult case of perturbations of the Kerr metric. As an example, we apply this method to the redshift invariant in Schwarzschild. Here we obtain analytic coefficients to 12.5PN, and higher-order terms in mixed analytic-numerical form to 21.5PN, including analytic forms for the complete 13.5PN coefficient, and all the logarithmic terms at 13PN. We have computed the individual modes to over 5000 digits, of which we use at most 1240 in the present calculation. At these high orders, an individual coefficient can have over 30 terms, including a wide variety of transcendental numbers, when written out in full. We are still able to obtain analytic forms for such coefficients from the numerical data through a careful study of the structure of the expansion. The structure we find also allows us to predict certain "leading logarithm"-type contributions to all orders. The additional terms in the expansion we obtain improve the accuracy of the PN series for the redshift observable, even in the very strongfield regime inside the innermost stable circular orbit, particularly when combined with exponential resummation.