A numerical method for solving the 2D TDGL equations in multiply connected domains has been described, together with the underlying term-by-term approximation formulae. An implementation of the method is available at http://cabmec1.cnea.gov.ar/gustavo, and is readily extendable to more than one hole to consider more complex mesoscopic systems. A few examples showing the main phenomena that can be simulated with the method have been reported. The reader is encouraged to download the program and test other configurations, such as eccentric holes (as considered in [15]) or squid-like systems. It is also possible to introduce defects, simply by changing in Eq. 3 by , with r<1 at regions with defects. Numerically, this amounts to a straightforward modification of Eq. 18. Extension to d-wave superconductors is more involved, but the term-by-term approximations of Section 2 lead to an appropriate scheme (some results can be found in [10,11]).