In this talk we present an algorithm for recursively generating orthogonal bivariate polynomials on a discrete set of the plane. For this purpose we employ commuting pairs of real symmetric matrices to obtain, in a certain sense, a two dimensional Hermitian Lanczos method. The resulting algorithm relies on a recurrence having a slowly growing length. Practical implementation issues and applications are considered. The method can be generalized to compute orthogonal polynomials depending on an arbitrary number of variables.