Statements on the behavior of solutions to $y' = \sin(xy)$ for large $y(0)$

Consider the following initial-value problem involving a nonlinear first-order ODE:

$y' = \sin(xy), \quad y(0) = y_0$.

For large enough values of $y_0$, the solutions to this equation appear to converge onto a specific shape under the appropriate normalized scale—dividing both $x$ and $y$ by $y_0$.

Can any statements be made about the form of the function corresponding to this shape/normalized solution?

See below for a graph of numerically computed solutions to this ODE that demonstrate this phenomenon.

I've plotted the solutions for even larger $y_0$ below so that the shape these solutions are converging to can be viewed more clearly. Only one line is visible as both solutions are virtually the same under the normalized scaling.

A comment: this question is motivated by a statement in Bender & Orszag's Advanced Mathematical Methods for Scientists and Engineers, where they use the IVP above when $y_0 = 1$ as an example of an initial-value problem with a unique solution near $x = 0$ that is not known to be representable in terms of elementary functions. Perhaps this strange shape is representable with elementary functions?

UPDATE 6/10/2018

Although it is not clearly visible in the graphs below, the function appears to have a slight kink at the point where $y = x$, indicating the solution's derivative is not itself differentiable. In addition, a paper by Wendell Mills, Boris Weisfeiler and Allan M. Krall discusses how solutions to $y' = sin(xy)$ asymptotically tend to a function $y \propto 1/x$ as $x \rightarrow \infty$.

Intrigued by this, I found numerically that for large $y_0$, the solution past the $y = x$ point (about 0.7929) matches $y = 0.63 y_0^2 / x$ essentially perfectly. See below for a comparison plot.

That being said, I still have no clue regarding:

1. What these constants—0.7929 and 0.63—are and where they come from. (I would assume the 0.63 constant has some factor of $\pi$ in it due to the arguments from the Mills-Weisfeiler-Krall paper.)
2. How to prove the solution past $x = y$ approaches this specific $\propto 1/x$ form for large $y_0$.
3. An expression for the part of the solution before $x = y$.

Even more intriguing: Perhaps the solution for large $y_0$ can be defined piecewise in terms of elementary functions?

UPDATE 7/31/2018

It seems as if the solution of $y' = \cos(xy),\ y(0) = y_0$ converges onto this unknown function as well for large $y_0$. See the graph below for numerical evidence; only one line is visible as both solutions are essentially the same.

Could the asymptotic form of these solutions have some exponential component in the region where $y > x$?

$$\frac{dy}{dx}=\sin(xy)$$

This is a very partial answer.

For large $x$ change of variable : $\quad x=\frac{1}{\epsilon}$ .

$$\frac{dy}{dx}=\sin(xy)=-\epsilon^2\frac{dy}{d\epsilon}=\sin(\frac{y}{\epsilon})$$ $y(\epsilon)=a_0+a_1\epsilon+a_2\epsilon^2+a_3\epsilon^3+a_4\epsilon^4+O(\epsilon^5)$ $$-\epsilon^2\left(a_1+2a_2\epsilon+3a_3\epsilon^2+4a_4\epsilon^3+O(\epsilon^4)\right)=\sin\left(\frac{a_0+a_1\epsilon+a_2\epsilon^2+a_3\epsilon^3+a_4\epsilon^4+O(\epsilon^5)}{\epsilon}\right)$$ $$-a_1\epsilon^2-2a_2\epsilon^3-3a_3\epsilon^4+O(\epsilon^5)=\sin\left(\frac{a_0}{\epsilon}+a_1+a_2\epsilon+a_3\epsilon^2+a_4\epsilon^3+O(\epsilon^4)\right)$$ This implies $\quad\begin{cases}a_0=0\\a_1=n\pi\end{cases}$ $$-a_1\epsilon^2-2a_2\epsilon^3-3a_3\epsilon^4+O(\epsilon^5)=(-1)^n\sin\left(a_2\epsilon+a_3\epsilon^2+a_4\epsilon^3+O(\epsilon^4)\right)$$ This implies $a_2=0$ . $$-a_1\epsilon^2-3a_3\epsilon^4+O(\epsilon^5)=(-1)^n\sin\left(a_3\epsilon^2+a_4\epsilon^3+O(\epsilon^4)\right)$$ $$-a_1\epsilon^2-3a_3\epsilon^4+O(\epsilon^5)=(-1)^n\left(a_3\epsilon^2+a_4\epsilon^3\right)+O(\epsilon^4)$$ This implies $-a_1=(-1)^na_3\quad$ and $\quad a_4=0$ . $$a_3=(-1)^{n+1}n\pi$$ $$y(\epsilon)=n\pi\epsilon+(-1)^{n+1}n\pi\epsilon^3+O(\epsilon^5)$$ $$y(x)=\frac{n\pi}{x}+(-1)^{n+1}\frac{n\pi}{x^3}+O\left(\frac{1}{x^5}\right)$$ This show that for $x\to\infty$ the function $y(x)$ tends to behave like $\frac{1}{x}$ , with a coefficient which is necessary an exact multiple of $\pi$.

With this method involving asymptotic series which cannot represent correctly the behavior for $x$ not large, there is no hope to determine $n$ as a function of $y_0$. This is a much more difficult part of the problem.

Note : I cannot open the paper by Wendell Mills, Boris Weisfeiler and Allan M. Krall which is probably more advanced than my small contribution.