The existence of crack and notch is a considerable and serious subject in the failure study of solids and structures. As most of cracked domain problems do not have closed-form solutions, numerical methods are dominant approaches dealing with fracture mechanics problems. In fracture mechanics and failure analysis, cracked media energy and consequently stress intensity factors (SIFs) play a crucial and significant role. Based on linear elastic fracture mechanics, the SIFs may be estimated using the energy release rate, as an indirect approach. This study presents the development of a new semi-analytical method to model cracked media. In this method, only the boundaries of problems are discretized using specific higher-order subparametric elements and higher-order Chebyshev mapping functions. Implementing the weighted residual method and using Clenshaw-Curtis quadrature yield diagonal Euler's differential equations. Therefore, when the local coordinate origin is located at the crack tip, the geometry of crack problems is implemented directly without further processing. In order to present crack boundary conditions, a modified form of nodal force function is presented. Validity and accuracy of the proposed method is fully illustrated by four benchmark problems, whose results agree very well with those of other numerical and/or analytical solutions existing in the literature.