We consider a boundary value problem (BVP) in unbounded 2D doubly periodic composite with circular inclusions having arbitrary constant conductivities. By introducing complex potentials, the BVP for the Laplace equation is transformed to a special -linear BVP for doubly periodic analytic functions. This problem is solved with use of the method of functional equations. The -linear BVP is transformed to a system of functional equations. A new improved algorithm for solution of the system is proposed. It allows one not only to compute the average property but to reconstruct the solution components (temperature and flux) at an arbitrary point of the composite. Several computational examples are discussed in details demonstrating high efficiency of the method. Indirect estimate of the algorithm accuracy has been also provided. Document embargo 02/01/2016.