The subject of the paper is the application of the Hilbert space version of the gradient method (GM) to linear elliptic boundary value problems of 2n-th order. In this setting the GM is applied to the generalized differential operator in the corresponding Sobolev space, thus reducing the original problem to auxiliary Poisson equations. The principle of realization is based on László Czách's method in the case of domains that can be transformed to a ball. Namely, the approximating sequence is constructed to consist of polynomials. This yields that the solution of the auxiliary equations can be achieved by solving systems of linear algebraic equations (SLAEs) with sparse matrices. The method yields linear convergence in Sobolev norm, moreover, for smooth enough coefficients we have uniform convergence. Numerical performance is investigated with the two-dimensional case in focus. The method is easy to realize, requiring the solution of SLAEs of simple structure, and yields linear convergence.