I didn’t get you when you say stddev is much larger than 1. The stddev I am using is 0.2 which translates to variance of 0.04. I agree that due to large stddev it may quickly decrease to zero or turn negative. So I tried with stddev of 0.05. It still does not work unless I use sqrt(abs(u)).
I replicated the problem in excel using sqrt(u) and it does not throw any negative there. The equation in excel is approximated as
r(t)=r(t-1)+a(b-r(t-1))*dt+sigma*sqrt(r(t-1))*sqrt(dt)*normsinv(rand())
I am not an expert in using DiffEq tool so may be missing something here.
Update
If I change sqrt(u) to sqrt(u*dt) in the equation even with sigma=0.2, it seems to work.