Let me add on the Feller condition.
With the sqrt model, there are 3 boundary regimes, which are the following (Feller scheme) classifications for the boundary of the V-process at V=0:
1. exit: [$] \omega \le 0[$]
2. regular: [$] 0 < \omega < \xi^2/2[$]
3. entrance: [$]\omega \ge \xi^2/2[$]
The posted code is valid for both regular and entrance regimes, while the Feller condition forces you to stay in the entrance regime. In the regular case, the interpretation of the solution is that it describes a V-particle which sometimes hits the origin and just reflects back to positive values. The posted code solution smoothly transitions from the unique solution in the entrance regime to its natural partner (reflecting) solution in the regular regime. There is lots more discussion in `A Closer Look at the Square-root and 3/2 model', which is Chapt . 7 in "Option Valuation under Stochastic Volatility II".