The problem of bivariate polynomial interpolation using Newton-type bases is examined, leading to the application of a generalized Kronecker matrix product. Algorithms for computing the coefficients of the interpolating polynomial are presented, along with conditions that ensure relative errors of the order of machine precision. A generalization of the classical recursion formula of divided differences in two dimensions is proposed for grids that generalize the standard rectangular layout. Numerical experiments demonstrate the high accuracy achieved by the proposed approach.