Magnetic induction tomography (MIT) attempts to image the passive electromagnetic properties (PEP) of an object by measuring the mutual inductances between pairs of coils placed around its periphery. In recent years, there has been an increase in applications of non-contact magnetic induction tomography. When finite element-based reconstruction methods are used, that rely on the inversion of a derivative operator, the large size of the Jacobian matrix poses a challenge since the explicit formulation and storage of the Jacobian matrix could be in general not feasible. This problem is aggravated further in applications for example when the number of coils is increased and in three-dimension. Krylov subspace methods such as conjugate gradient (CG) methods are suitable for such large scale inverse problems. However, these methods require use of the Jacobian matrix, which can be large scale. This paper presents a matrix-free reconstruction method, that addresses the problems of large scale inversion and reduces the computational cost and memory requirements for the reconstruction. The idea behind the matrix-free method is that information about the Jacobian matrix could be available through matrix times vector products so that the creation and storage of big matrices can be avoided. Furthermore the matrix vector multiplications were performed in multiple core fashion so that the computational time can decrease even further. The method was tested for the simulated and experimental data from lab experiments, and substantial benefits in computational times and memory requirements have been observed.