Update to follow-up

Let me rewrite your model one more time $$ Y_i = \beta_0 + \beta_1Z_{1i} + \beta_2 Z_{2i} + \beta_3 X_i + U_{1i} $$ where as before, $Z_{1i}$ and $Z_{2i}$ are exogenous and $X_i$ is endogenous. You assume the data generating process $$ \begin{align} Z_{2i} &= \mathbf{1}_{[i\text{ is odd}]}\\ X_i &= \alpha_0 + \alpha_1Z_{1i} + U_{2i} \end{align} $$ You introduce endogeneity by assuming that $U_{1i}$ and $U_{2i}$ are correlated. $$ \mathbb{C}(U_{1i}, U_{2i}) = \rho $$

Also, you assume that the errors in the reduced form equation are heteroskedastic: $$ \begin{align} \mathbb{V}(U_{2i}\mid Z_{2i} = 0) &= 1 \\ \mathbb{V}(U_{2i}\mid Z_{2i} = 1) &= \tfrac{1}{q} \\ \end{align} $$

Here is the Stata code to simulate this DGP.

clear*
program simcont, rclass
    syntax [, Q(real 1.0)]

    drop _all
    set obs 1000
    scalar beta0 = 5
    scalar beta1 = 1
    scalar beta2 = 2
    scalar beta3 = 3

    scalar alpha0 = 1
    scalar alpha1 = 4

    scalar rho = 0.1

    g z1 = rnormal()
    g z2 = mod(_n, 2)
    scalar a11 = 1
    scalar a12 = rho*sqrt(1)*sqrt(1)
    scalar a13 = rho*sqrt(1)*sqrt(1/`q')
    scalar a21 = rho*sqrt(1)*sqrt(1)
    scalar a22 = 1
    scalar a23 = 0
    scalar a31 = rho*sqrt(1)*sqrt(1/`q')
    scalar a32 = 0
    scalar a33 = 1/`q'

    mat corrMatrix= (a11, a12, a13 \ a21, a22, a23 \a31, a32, a33)
    drawnorm u1 u3 u4, cov(corrMatrix)
    g u2 = cond(z2, u3, u4)

    g x = alpha0 + alpha1*z2 + u2
    g y = beta0 + beta1*z1 + beta2*z2 + beta3*x + u1
    reg y z1 z2 x
    mat mA = e(b)
    return scalar ols = el(mA, 1, colnumb(mA, "x"))
    
    // return the results of the heteroskedasticity test
    qui reg x z2
    estat hettest, rhs iid
    return scalar hettestPValues = r(p)
end

// simulate to check for heteroskedasticity
simulate olsCoeff = r(ols) hettestPValues = r(hettestPValues), reps(100): simcont, q(5) 
su hettestPValues // strong evidence of heteroskedasticity in the reduced form

cap mat drop biasBeta
forvalues q = 2(1)5 {
    simulate olsCoeff = r(ols), reps(1000): simcont, q(`q') 
    qui su olsCoeff
    mat biasBeta = (nullmat(biasBeta), r(mean) - 3)
    local colNames "`colNames' q:`q' "
}
mat colnames biasBeta = `colNames'
mat list biasBeta

• Note that the heteroskedasticity is where you put it in the model, and the simulation results are now able to find it.

. su hettestPValues // strong evidence of heteroskedasticity in the reduced form

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
hettestPVa~s |       100    1.53e-26    8.12e-26   3.09e-38   6.22e-25

• Also, I can confirm that what you claim is actually true, that as the heteroskedasticity in the reduced form equation increases, the bias in the OLS estimator also increases.

biasBeta[1,4]
            q:         q:         q:         q:
            2          3          4          5
r1  .11242463  .11808594  .12122701  .12210455

A simple explanation of this is that as $q$ increases, the conditional (on $Z_{2i}=1$) and hence the unconditional variance of the reduced form error $U_{2i}$ decreases (check through a variance decomposition) which means that the variance of the endogenous regressor decreases, which means that the $(\mathbf{X}'\mathbf{X})^{-1}$ increases in magnitude, and the overall bias increases (this is a very rough description -- I am sure it can be formalized).

Follow-up questions:

The OP has asked for some clarifications in the comments for which I am providing an updated answer.

Let me rewrite your model $$ Y_i = \beta_0 + \beta_1 Z_{1i} + \beta_2 Z_{2i} + \beta_3 X_i + U_i $$ where $Z_{1i}$ and $Z_{2i}$ are exogenous, and $X_{i}$ is endogenous, that is $\mathbb{E}(X_iU_i) = \rho$. In addition, you want the variable $Z_{2i}$ to be constructed as $$ Z_{2i} = \mathbf{1}_{[i\text{ is odd.}]} $$ You are simulating the OLS estimates of the coefficient on $X_i$, that is $\beta_3$. Here is a small Stata script to do that.

clear*
program simcont, rclass
    syntax [, errorVariance(real 1.0)]
    drop _all
    set obs 1000
    scalar beta0 = 5
    scalar beta1 = 1
    scalar beta2 = 2
    scalar beta3 = 3
    scalar rho = 0.1

    g z1 = rnormal()
    g z2 = mod(_n, 2)
    g u = sqrt(`errorVariance')*rnormal() 
    g x = rho*u/`errorVariance' + rnormal()

    g y = beta0 + beta1*z1 + beta2*z2 + beta3*x + u
    reg y z1 z2 x
    mat mA = e(b)
    return scalar ols = el(mA, 1, colnumb(mA, "x"))
    
    // return the results of the heteroskedasticity test
    estat hettest, rhs iid
    return scalar hettestPValues = r(p)
end

// simulate with error variance = 1
simulate olsCoeff = r(ols) hettestPValues = r(hettestPValues), reps(100): simcont 
su olsCoeff hettestPValues // p-values have the correct mean; no heteroskedasticity

cap mat drop biasBeta
forvalues errorVariance = 2(1)6 {
    simulate olsCoeff = r(ols), reps(100): simcont, errorVariance(`errorVariance') 
    qui su olsCoeff
    mat biasBeta = (nullmat(biasBeta), r(mean) - 3)
    local colNames "`colNames' errVar:`errorVariance' "
}
mat colnames biasBeta = `colNames'
mat list biasBeta

Now note that there are at least two discrepancies here.

• The first thing is that there is no heteroskedasticity in the model as you have written it. The mean of the p-values of an LM test of homoskedasticity from the simulations indicates that they are being drawn under the null.

. su olsCoeff hettestPValues // p-values have the correct mean; no heteroskedasticity

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
    olsCoeff |       100    3.094182    .0326048   3.007202   3.179945
hettestPVa~s |       100    .4762952    .2885692   .0039824   .9844939

• The next thing to note is that the bias is that the asymptotic bias is the same for all values of the error variance, as long as the degree of endogeneity is the same. Here are the results of my simulation

. mat list biasBeta

biasBeta[1,5]
       errVar:    errVar:    errVar:    errVar:    errVar:
            2          3          4          5          6
r1  .09913641  .10274499   .0912232  .11309942  .09122764

If you are getting results different than this, then you should show us your code, and we can compare the two.