If the nonhomogeneous part f of the linear elliptic PDE in the region D with boundary C, Lu= f in D, Bu= g on C, can be approximated by a polynomial in its variables x i , i=1,…, n, n=2, (3), then it is possible to write a particular solution p u of this equation as a linear combination of the particular solutions P kl(m) of Lu= x k 1 x l 2( x m 3). The double (triple) sequence of solutions P kl(m) can be determined rec consequence. The particular solution u p of the given equation can also be evaluated recursively in each point of interest. This allows us to reduce the given nonhomogeneous problem to the homogeneous case, and to use the recursive method already described for that problem to find a solution u h of this new elliptic problem, Lu h=0 in D, Bu h= g- Bu p on C, and thus obtaining the solution u of the initial problem as a sum of u p and u h. It is to be emphasized that this solution, as well as its derivatives, can be evaluated in each point from the coefficients of u p and u h in a fast and stable way.