In this paper, a numerical method is developed to compute an approximate solution of high-order linear complex differential equations in elliptic domains. By using collocation points and Bessel polynomials, this method transforms the linear complex differential equation into a matrix equation which corresponds to a system of linear equations in the unknown Bessel coefficients. The suggested method reproduces the analytic solution when it is polynomial. Some numerical examples are given to illustrate the applicability and validity of the technique and comparisons are made with other existing approaches. © de Gruyter 2011.