The Fortran implementation of an iterative lookahead algorithm for the computation of a diagonal of the Padé tables of Padé-Hermite and simultaneous Padé systems is described. The results of the computation include certain Padé-Hermite and simultaneous Padé approximants that are defined as vectors of polynomials subject to some order and nonsingularity conditions. These approximants can be regarded as generalizations of the ordinary Padé approximant for a single power series to vectors of power series. The computation requires the solution of some system of linear equations of order O ( N ), where N is the sum of the (maximal) degrees of the polynomials of the approximant. The algorithm is O ( N ² ) arithmetically, that is, faster than Gaussian elimination due to a structured form of the coefficient matrix, and weakly stable, in contrast to other fast and superfast O ( N ln² ( N ) ) algorithms for this problem. This weak stability is achieved by avoiding badly conditioned elements in the Padé tables, by using a stability parameter based on the nonsingularity conditions. This lookahead criterion can make the algorithm O ( N⁴ ) in the worst case. As an illustration, a basic step of the iterative algorithm is described and a worked example is given. Some numerical evidence shows that the algorithm works and that several error bounds are satisfied, albeit too pessimistically. This rather long research paper is useful for experts designing rational approximation algorithms, but not for casual readers, because of the technical presentation and the lack of practical applications.