Let's make a substitution $\sin(\theta)=y$
$$
I(x,z)=\int_{0}^{\sin(z)}\frac{J^2_{1}(x y)}{y^2\sqrt{1-y^2}}dy
$$
Because $z\approx0.1$ we may that assume $\sin(z)\approx z$. Furthermore we perform another subsitution $\frac{y}{z}=q$ yielding
$$
I(x,z)=z\int_{0}^{1}\frac{J^2_{1}(x z q)}{(zq)^2\sqrt{1-(zq)^2}}dq
$$
As a first step to simplify the problem we may observe that for $z\approx 0.1$, $\sqrt{1-(zq)^2}\approx 1$ over the whole domain of integration. Therefore:
$$
I(x,z)\approx\frac{1}{z}\underbrace{\int_{0}^{1}\frac{J^2_{1}(x z q)}{q^2}dq }_{Q(x,z)}\quad(1)
$$
Essentially there are two ways to go from here on:
1.) Ask Mathematica what the result of this integral is:
$$
I(x,z)\approx\frac{1}{4} x^2 z \, _1F_2\left(\frac{1}{2};2,3;-(x z)^2\right) \quad(2)
$$
Note that this result is valid for all values of $x$! only the assumption that $z$ is sufficiently smaller then 1 was made!. If you are interested in the large $xz$ you may use 15.8.2 from here which connects small arguments of $_2F_1$ to big arguments of $_2F_1$. If i have done everything correctly this matches with the answer found by @Raymond Manzoni
$$
I(x,z)\sim\frac{4}{3\pi}x
$$
2.) Second approach would be to really to dig in (1) and derive the explicit results for $xz \rightarrow \infty$ using a split of integration range and asymptotic expansions of Bessel functions. I hope i can finish this approach today or tomorrow when i'm less busy. But i'm quite sure this works out fine.
$\bf{Appendix}$
To prove (2) we may employ the following representation of $J^2_1(x)$:
$$
J^2_1(x)=\sum_{m=0}^{\infty}\frac{(-1)^m}{2^{2+2m}}\frac{x^{2m+2}\Gamma(2m+3)}{m!\Gamma(m+3)\Gamma^2(m+2)}
$$
Now after switching summation and integration
$$
Q(x,z)=\sum_{m=0}^{\infty}\frac{(-1)^m}{2^{2+2m}(2m+1)}\frac{(zx)^{2m+1}\Gamma(2m+3)}{m!\Gamma(m+3)\Gamma^2(m+2)}
$$
By application of the duplication formula of the gamma function this may be reduced to
$$
Q(x,z)= \sum_{m=0}^{\infty}\frac{(-1)^m}{2\sqrt{\pi}}\frac{(zx)^{2m+1}\Gamma(m+\frac{1}{2})}{m!\Gamma(m+3)\Gamma(m+2)}
$$
Using the definition of Pochhammer symbol $(x)_n=\frac{\Gamma(x+n)}{\Gamma(x)}$ we find that $(\frac{1}{2})_n=\sqrt{\pi}\Gamma(n+\frac{1}{2}),(2)_n=\Gamma(n+2)$ and $(3)_n=2\Gamma(n+2)$ so we can reformulate
$$
Q(x,z)= \frac{z^2x^2}{4}\sum_{m=0}^{\infty} \frac{(-z^2x^2)^{m}(\frac{1}{2})_m}{m!(2)_m (3)_m}
$$
Using the series representation of $_1F_2$ we may conclude that
$$
Q(x,z)=\frac{z^2 x^2}{4}{_1F}_2\left(\frac{1}{2};2,3;-(x z)^2\right)
$$
Which proves (2)