We provide a comparative analysis of qualitative features of different numerical methods for the inhomogeneous geometric Brownian motion (IGBM). The conditional and asymptotic mean and variance of the IGBM are known and the process can be characterised according to Fellers boundary classification. We compare the frequently used Euler-Maruyama and Milstein methods, two Lie-Trotter and two Strang splitting schemes and two methods based on the ordinary differential equation (ODE) approach, namely the classical Wong-Zakai approximation and the recently proposed log-ODE scheme. First, we prove that, in contrast to the Euler-Maruyama and Milstein schemes, the splitting and ODE schemes preserve the boundary properties of the process, independently of the choice of the time discretisation step. Second, we derive closed-form expressions for the conditional and asymptotic means and variances of all considered schemes and analyse the resulting biases. While the Euler-Maruyama and Milstein schemes are the only methods which may have an asymptotically unbiased mean, the splitting and ODE schemes perform better in terms of variance preservation. The Strang schemes outperform the Lie-Trotter splittings, and the log-ODE scheme the classical ODE method. The mean and variance biases of the log-ODE scheme are very small for many relevant parameter settings. However, in some situations the two derived Strang splittings may be a better alternative, one of them requiring considerably less computational effort than the log-ODE method. The proposed analysis may be carried out in a similar fashion on other numerical methods and stochastic differential equations with comparable features.