How can the increments of a CIR process be derived?

I am not totally sure I understand what you want to achieve. It seems like you are interested in discretizing CIR SDE. This can be done using the Euler-Murayama scheme for an equidistant decomposition of the time interval $[0, T]$, ${0=t_0<dots<t_n=T}$.

First of all, let us write the model dynamics: $$r_t=r_0+alphaint_0^t(mu-r_s)ds+sigmaint_0^tsqrt{r_s}dW_s$$

We need to discretize this process: $$r_{t+Delta t}=r_t+alpha(mu-r_t)Delta t+sigmasqrt{r_t}W_{Delta t}$$ with $Delta t=frac{T}{n}$ and $W_{Delta t}simmathcal Nleft(0,frac{T}{n}right)Rightarrow W_{Delta t}=sqrt{frac{T}{n}}varepsilon,$ with $varepsilon$ being a standard normal random variable.

Finally, we can use the trapezoidal rule to numerically integrate the simulated CIR rates and compute what you need (for example, the Monte Carlo zero-coupon bond prices).