Exact comprehensive computations are carried out by means of four leading second-order approximations yielding differential cross sections dQ/d Omega for the basic charge exchange process H+ + H(1s) --> H(1s) + H+ at intermediate and high energies. The obtained extensive set of results is thoroughly tested against all the existing experimental data with the purpose of critically assessing the validity of the boundary corrected second-Born (CB2), continuum-distorted wave (CDW), impulse approximation (IA) and the reformulated impulse approximation (RIA). The conclusion which emerges from this comparative study clearly indicates that the RIA agrees most favorably with the measurements available over a large energy range 25 keV-5 MeV. Such a finding reaffirms the few-particle quantum scattering theory which imposes several strict conditions on adequate second-order methods. These requirements satisfied by the RIA are: (i) normalisations of all the scattering wave functions, (ii) correct boundary conditions in both entrance and exit channels, (iii) introduction of a mathematically justified two-center continuum state for the sum of an attractive and a repulsive Coulomb potential with the same interaction strength, (iv) inclusion of the multiple scattering effects neglected in the IA, (v) a proper description of the Thomas double scattering in good agreement with the experiments and without any unobserved peak splittings. Nevertheless, the performed comparative analysis of the above four approximations indicates that none of the methods is free from some basic shortcomings. Despite its success, the RIA remains essentially a high-energy model like the other three methods under study. More importantly, their perturbative character leaves virtually no room for further systematic Improvements, since the neglected higher-order terms are prohibitively tedious for practical purposes and have never been computed exactly. To bridge this gap, we presently introduce the variational Fade approximant (VPA) as a novel non-perturbative theory which is valid at all energies. This is a variationally unified T-matrix, T-i,f((VPA)) = T-i,T-f' + S-i,S-f comprised of a selected perturbative model T-i,T-f' and a related stationary non-linear remainder S-i,S-f conceived as an L-2-basis set expansion of the total Green's function. The key input S-i,S-f is a double series Sigma(n1,n2) C-n1n2 With bound-free atomic form factors as rational coefficients C-n1n2. Convergence of this sum is significantly accelerated by the two-dimensional Pade approximant OD-PA) implemented through the bi-variate Wynn's epsilon (epsilon) algorithm. The table epsilon(lambda) for S-i,S-f(lambda) is evaluated at a sufficiently dense grid lambda is an element of Lambda = [lambda(min), lambda(max)] of any chosen variate lambda, e.g. scattering angle, incident energy, coupling strength, etc. The roots lambda(k) of the inverse function E(lambda) = 1/epsilon(lambda) directly lead to the poles of S-i,S-f(lambda) in the complex lambda-plane. The availability of all lambda(k)'s permits an easy obtaining of each of the magnitudes d(k) = 1/E'(lambda(k)) of the poles of S-i,S-f(lambda) by the Cauchy residue calculus and a mere knowledge of the first derivative E'(lambda(k)) of E(lambda(k)) at lambda(k) is an element of Lambda. This is a new signal processing method called the parametric epsilon spectral (PES) estimator which opens an affordable way to powerful numerical investigations of analytical properties of transition and scattering matrices for a general process including resonance phenomena as one of the the most interesting parts of scattering. (C) 1999 Elsevier Science B.V. All rights reserved.