The modified Bessel differential equation can be obtained by replacing $x$ with $ix$ ($i$ is the imaginary unit) in the Bessel differential equation. If the general solution of the latter is
$$f(x) = A J_{n}(x) + BN_n(x)$$
While the general solution of the modified Bessel differential equation is
$$g(x) = CI_n(x) + DK_n(x)$$
Knowing the definitions of $I_n(x)$ and $K_n(x)$, I don't find a way to express $g(x)$ as $f(ix)$. That is, as far as I know,
$$f(ix) = A J_{n}(ix) + BN_n(ix) \neq g(x) = CI_n(x) + DK_n(x)$$
But if the modified Bessel differential equation is simply obtained by substituting the variable $x$ with $ix$ in the Bessel differential equation, why this can not be applied to the general solutions too?
In other words, why the change of variable can not immediately be "transferred" also to the solutions?
The behavior of solutions of $$ r^2f''(r)+rf'(r)+(r^2-\nu^2)f(r) = 0 $$ near $r=0$ is asymptotically the same as the solutions of $$ r^2f''+rf'-\nu^2 f=0. $$ For $\nu > 0$, the solutions of the latter are $r^{\pm\nu}$. For $\nu=0$, the solutions are $1$ and $\ln(r)$. Because of this, solutions in the Complex plane inherit branch points at the origin, and the conventions used to deal with branch cuts are not uniquely determined. Solutions are locally holomorphic, but somewhere there must be a branch cut if you work in the Complex plane.
Example of Asymptotic for $\nu=0$: In this case, the solutions of $$ r^2f''+rf'=0 $$ are $1$ and $\ln r$. You can check that $$ -\frac{1}{r^2}\left[r^2\frac{d^2}{dr^2}+r\frac{d}{dr}\right] \left(A+B\ln(r)+\int_{0}^{r}g(s)\ln(s)sds-\ln(r)\int_{0}^{r}g(s)sds\right)=g $$ (I may have the signs in front of the two integral terms swapped, but I'll let you work that out. It's a Green function solution.) Therefore a solution of $$ \left[r^2\frac{d^2}{dr^2}+r\frac{d}{dr}\right]f+r^2f=0 \\ -\frac{1}{r^2}\left[r^2\frac{d^2}{dr^2}+r\frac{d}{dr}\right]f=f $$ must satisfy the following for some constants $A$ and $B$: $$ A+B\ln(r)+\int_{0}^{r}f(s)\ln(s)sds-\ln(s)\int_{0}^{r}f(s)sds=f(r) $$ The integral terms must vanish at $r=0$ because the equation is in the limit circle case at $r=0$. This gives an asymptotic of $f\sim A+B\ln(r)$ near $r=0$. This is the simplest asymptotic, and this technique works for $0 \le \nu < 1$ because the equation is in the limit point case at $r=0$ for all such $\nu$, meaning that all eigenfunctions are in $L^2_{r}(0,1]$ (the weighted space $(f,g)_r=\int_{0}^{1}f(r)\overline{g(r)}rdr$.)