On the smooth dependence on initial conditions for the following ODE

The continuous dependence of the solution on the initial value is proven using the

Grönwall inequality: for $x>x_0$ one has $$ \|\Delta u'(x)\|\le L(x)\,\|Δu(x)\|\implies \|Δu(x)\|\le\exp\left(\int_{x_0}^xL(s)ds\right)\,\|Δu(x_0)\|. $$

For this problem use $u=(y,y')$, $\|u(x)\|=\max(|y(x)|,|y'(x)|)$ so that then on $x>0$ $$ \|\Delta u'(x)\|\le\max(|Δy'(x)|,\frac3{2x}|Δy(x)|+x^{7/2}|Δy(x)|)\le\max(1,\frac3{2x}+x^{7/2})\,\|Δu(x)\|. $$

So essentially, there are 2 segments for the Lipschitz constant,

• for $x<1.1$, the Lipschitz constant is dominated by $L(x)\simeq \frac3{2x}$ leading to $$\|Δu(x)\|\le\left(\frac{x}{x_0}\right)^{3/2}\,\|Δu(x_0)\|$$
• for $x>1.1$ the Lipschitz constant is dominated by the second term $L(x)\simeq x^{7/2}$ so that $$ \|Δu(x)\|\le\exp\left(\frac29(x^{9/2}-1.1^{9/2}\right)\|Δu(1.1)\| \le\exp\left(\frac29(x^{9/2}-1.1^{9/2}\right)\left(\frac{1.1}{x_0}\right)^{3/2}\,\|Δu(x_0)\| $$

As one can see, on the second segment the factor relating the error at $x_0$ to the error at $x$ deteriorates rapidly, for instance to get $\|Δu(5)\|<1$ one would need at $x_0=0.1$ $$ \|Δu(0.1)\|\le \exp\left(-\frac29(5^{9/2}-1.1^{9/2}\right)\left(\frac{0.1}{1.1}\right)^{3/2}\le 5.12205\cdot 10^{-137} $$ which is impossible to realize in standard floating point values. Even allowing that this estimate is somewhat pessimistic, the practical result will be a jump in behavior for minimal changes in the initial conditions.

---

To make an observable prediction, let the initial condition have a variance in the machine epsilon range, $\|Δu(x_0)\|\sim 10^{-15}$ and calculate the interval where the solutions are visually close, that is $\|Δu(x)\|<10^{-2}$. There is no problem on the first segment, so on the second segment we have to solve $$ \exp\left(\frac29(x^{9/2}-1.1^{9/2}\right)\left(\frac{1.1}{0.1}\right)^{3/2} \le 10^{13}$$ leading to $x\le2.9$. This seems consistent with the plots in the question as the graphs start to be different at about $x=3.5$.

def sol(x,y0): return odeint(lambda y,x:[y[1],-3/(2*x)*y[0]-x**3.5*np.sin(y[0])], y0, x, atol=1e-12, rtol=1e-13) 

x = np.linspace(0.1,5,703);
for k in range(-4,5):
    a = 2.2191851423+1e-7*k; 
    y=sol(x,[a, 1]); 
    plt.plot(x,y[:,0]);
plt.grid(); plt.show()

As one can see, the actual behavior is several orders of magnitude better than predicted, the non-linearity appears to also work towards stabilizing the solution behavior towards shrinking oscillations around multiples of $2\pi$, so that one does not see the worst case estimate in every interval of initial values.