Inclusions in engineering materials will create a disturbance to the elastic fields. When an inclusion is located in the vicinity of the surface, the interactions between the inclusion and the surface are often difficult to be solved analytically. The elementary solutions of the displacements and stresses caused by a rectangular inclusion contained in an elastic half-plane are derived in closed-form. Accordingly, the resultant elastic fields produced by any arbitrarily shaped inclusion may be solved through a "discretization-superposition" scheme. In contrast, the traditional Finite Element Method needs to mesh the semi-infinite matrix region which is much larger than the size of the inclusion, while the mesh needs to be sufficiently refined at the inclusion/matrix interface to attain a satisfactory accuracy. By taking advantages of the elementary solutions of rectangular inclusion, the proposed semi-analytical method merely performs the numerical discretization within the inclusion region, leading to remarkable savings on the meshing efforts and memory storage. Benchmark examples are reported for an elastic half-plane containing either regular hexagonal or circular inclusion, and the results are validated with those obtained by the commercial Finite Element software to show the correctness and effectiveness of this proposed semi-analytical method.