Summary: We describe a new method of calculation of generic multiloop master integrals based on the numerical solution of systems of difference equations in one variable. We show algorithms for the construction of the systems using integration-by-parts identities and methods of solutions by means of expansions in factorial series and Laplace transformation. We also describe new algorithms for the identification of master integrals and the reduction of generic Feynman integrals to master integrals, and procedures for generating and solving systems of differential equations in masses and momenta for master integrals. We apply our method to the calculation of the master integrals of massive vacuum and self-energy diagrams up to three loops and of massive vertex and box diagrams up to two loops. Implementation in a computer program of our approach is described. Important features of the implementation are: the ability to deal with hundreds of master integrals and the ability to obtain very high precision results expanded at will in the number of dimensions.