In this thesis a new boundary point method is presented for the solution of problems in solid mechanics. The method is based on the use of basis functions satisfying the governing partial differential equations (PDEs) approximately. The idea is in contrast with the approach usually employed in the conventional methods using fundamental solutions. In this method the basis functions are found through satisfaction of the weighted residual of the PDEs written in integral forms. The primary bases and the weight functions used for obtaining the final basis functions are, respectively, considered as Chebyshev polynomials (of the first kind) and exponential functions. The reason of this selection is the consistency between the basis and the weight functions which gives good characteristics to the method. The integrals are taken on a fictitious rectangular domain circumscribing the actual domain. This helps to compute the multidimensional integrals by multiplication of simpler one-dimensional ones. This is possible because the basis and the weight functions can always be written as multiplication of a function of x into a function of y. Such decomposition also speeds up the computational procedure significantly. After forming all the required integrals, they will be arranged in a matrix equation to compute the constant coefficients of solution series. By solving the matrix equation, a set of secondary basis functions will be created that satisfy the homogeneous operator approximately. The new bases form another solution series with new unknown coefficients which should be obtained by applying the boundary conditions. A special discrete transformation is used for this purpose which calculates each coefficient according to the projection of the corresponding basis on the boundary. Different shapes of domain can be considered in the method. The method is capable of solving PDEs with non-constant coefficients. A special technique is developed for this part so as to continue the use of separation method discussed for PDEs with constant coefficients. In this technique the so called Pascal triangle is used to rewrite the non-constant coefficients of operators as the sum of simple polynomials. Then the one-dimensional integrals are calculated for each part of the series and combined through a special matrix operation. By using this technique, the speed of calculation for equations with constant and non-constant coefficient will be of the same order. Three partial differential operators are discussed in this research; namely operators in engineering problems known as: Helmholtz, elasto-static and elastic wave and Kirchhoff plate. Each operator is formulated in two approaches for the integral equation known as the strong form and the weak form. The latter has the advantage of removing all differential operators from the non-constant coefficients of the equation which may have complicated behavior. The formulation of thesis is for two-dimensional problems, but it is possible to extend it to three-dimensional cases while holding all advantages addressed for two-dimensional problems. The examples solved throughout the text confirm the capability and accuracy of the method in solution of PDEs with both constant and non-constant coefficients. Convergence study is given for all examples solved. Key Words Partial differential equation, Weighted residual integral, non-constant coefficients, Chebyshev polynomials