A meshless method based on the local Petrov-Galerkin approach is applied for boundary value problems with cracks in magnetoelectroelastic solids. A unit step function is used as the test functions in the local weak-form. This leads to local boundary-domain integral equations (LIEs). The moving least-squares (MLS) method is adopted for approximating the physical quantities in the LIEs. A system of ordinary differential equations for certain nodal unknowns is obtained. That system is solved numerically by the Houbolt finite-difference scheme. Numerical results for intensity factors in homogeneous and functionally graded materials are presented.