I have a diffusion process $X=(X_t)_{t \ge 0}$ with the generator $$Af(x)=\frac{1}{2}(a(1-x)-bx)f'(x)+\frac{1}{4}x(1-x)f''(x),$$ where $a,b >0$ are constants. I want to simulate $X$ to a deterministic time $T$ if I start from $X_0=x_0 \in [0,1]$. So the obvious way is to discretize it to equidistant times $X_{t_i},\,i=1,\ldots,N$, so that $t_i-t_{i-1}=\Delta := T/N$. Then we can approximate $X_{t_i}$ by $$X_{t_i}=X_{t_{i-1}}+\frac{1}{2}\Delta(a(1-X_{t_{i-1}})-bX_{t_{i-1}})+\frac{1}{\sqrt{2}}\sqrt{\Delta}Z_i\sqrt{X_{t_{i-1}}(1-X_{t_{i-1}})},$$
where $Z_i,\,i=1,\ldots,N$ are i.i.d. such that $Z_i \sim N(0,1).$ However, when I actually do that, then the simulation seems to depend on the choice of $\Delta$ and eventually exits the interval $[0,1]$ (which the process never leaves and actually since $a$ and $b$ are positive it never becomes trapped on the boundary either). This seems to be caused by the fact that the diffusion part is multiplied by $\sqrt{\Delta}$ whereas the drift only by $\Delta$ and so the diffusion part is much stronger. However, since for a Wiener process $W_t-W_s \sim N(0,t-s),$ this should to be correct.
I guess I am probably missing something trivial or else misunderstand something fundamental, but I cannot see an error in my reasoning while it obviously cannot be correct.