Random samples from a given 2D distribution

I'd like to generate random samples based on given 2D PMF $p=(\theta,\phi)$ with dependent variables.

For a 1D PMF, inverse transform sampling was done and verified simply by doing numerical integrations to compute CDF and then a simple interpolation (interp1 from matlab) for doing the inverse.

For 2D PMF, using this discussion, I think I need to compute conditional CDF by dividing marginal CDF (with respect to first variable) dividing by a joint CDF as mentioned

Then, using the marginal of X we can get the conditional distribution of Y|X by dividing the joint by the marginal Then integrate that with respect to y to get the conditional CDF of Y|X.

1. can the conditional CDF be computed by cdfCond=cdfJoint./cdfMarginal where cdfJoit and cdfMarginal are matrices, computed from 2d numerical integrations? I mean, is it a simple element-wise division? (The code for computing jointCDF can be found here)

2. The conditional and marginal CDF have to be inverted finally. A 1D (non-decreasing) function can be easily inverted using an interpolation function like "interp1" in matlab. But the problem is that my CDFs are now matrices. I am confused here, because it seems that i have to invert a 1D function as written here

Now to get a random sample from a uniform[0,1] for X & transform it

Xs<-runif(1000000)

Xs<-margXcdf_inv(Xs)