Composite structures are distinguished by interfaces separating their different layers. These interfaces represent region of mismatches in mechanical and thermal properties of the adjacent layers. These regions are rich in cracks which are laying either on or sub- the interface in a direction parallel or perpendicular to it. In bi-material composites, interfacial crack is under mixed loading mode and stress intensity factor is a complex number even if the opening mode loading is applied. This paper solves for energy release rates and the magnitude of complex intensity factor of interfacial cracks using the classical stiffness finite element method. The path-independent integral is used to calculate the energy release rate during crack extension in the elastic plane strain deformation field. In this work, the crack tip region is meshed with fan-shaped non-singular four nodes isoparametric quadrilateral elements. In addition, a conical shape of function is employed in the numerical implementation of the equivalent energy domain integral. Realization of all numerical computations and demonstration of results are completely composed and written in MATLAB language. Meshing the computational domains and crack tip region are performed by a free downloadable program AUTOMESH-2D. Numerical results of stress intensity factor are found to be very close to the analytical and referenced values in both cases of bi-material and single layer systems. Furthermore, numerical values of -integral contours are very close to interfacial fracture energy measured experimentally between a hard film and a soft substrate.