Thank You for your fast response, but why should an iterative method converge if the matrix is completely random? it's not positive, not diagonal dominant and not symmetric... Does your code work always?

Not always converge. In most of my cases, due to central difference used for Laplace operator, diagonal dominat is available.

Maybe you use the elliptic grid generator as smoothing for a pre existing grid; in this case you should find coefficients for the elliptic equation that keep the dominance.

But could you tell mer something about the orthogonality?