Polynomial coding has been proposed as a solution to the straggler mitigation problem in distributed matrix multiplication. Previous works employ univariate polynomials to encode matrix partitions. Such schemes greatly improve the speed of distributed computing systems by making the task completion time to depend only on the fastest workers. However, they completely ignore the work done by the slowest workers resulting in inefficient use of computing resources. In order to exploit the partial computations of the slower workers, we further decompose the overall matrix multiplication task into even smaller subtasks, and we propose bivariate polynomial codes. We show that these codes are a more natural choice to accommodate the additional decomposition of subtasks, and to exploit the heterogeneous storage and computation resources at workers. However, in contrast to univariate polynomial decoding, guarantying decodability with multivariate interpolation is much harder. We propose two bivariate polynomial coding schemes and study their decodability conditions. Our numerical results show that bivariate polynomial coding considerably reduces the computation time of distributed matrix multiplication.