We develop a new simulation method for multidimensional diffusions with sticky boundaries. The challenge comes from simulating the sticky boundary behavior, for which standard methods like the Euler scheme fail. We approximate the sticky diffusion process by a multidimensional continuous time Markov chain (CTMC), for which we can simulate easily. We develop two ways of constructing the CTMC: approximating the infinitesimal generator of the sticky diffusion by finite difference using standard coordinate directions, and matching the local moments using the drift and the eigenvectors of the covariance matrix as transition directions. The first approach does not always guarantee a valid Markov chain whereas the second one can. We show that both construction methods yield a first order simulation scheme, which can capture the sticky behavior and it is free from the curse of dimensionality. We apply our method to two applications: a multidimensional Brownian motion with all dimensions sticky which arises as the limit of a queuing system with exceptional service policy, and a multi-factor short rate model for low interest rate environment in which the stochastic factors are unbounded but the short rate is sticky at zero.