A recently developed numerical method based on a volume integral formulation is applied to calculate the accurate stress intensity factors at the crack tips in unbounded isotropic solids in the presence of multiple anisotropic inclusions and cracks subject to external loads. In this paper, a detailed analysis of the stress intensity factors are carried out for an unbounded isotropic matrix containing an orthotropic cylindrical inclusion and a crack. The accuracy and effectiveness of the new method are examined through comparison with results obtained from analytical method and finite element method using ANSYS. It is demonstrated that this new method is very accurate and effective for solving plane elastostatic problems in unbounded solids containing anisotropic inclusions and cracks.