Yes, but.
The quality of the collimated beam may not be fantastic. And it depends on how the light is launched into the fiber. If the incident light has a narrow angular dispersion, and the fiber is short and/or of very good quality with few bends, then the light coming out will have a low angular dispersion. But if the incident light is converging at an angle close to the NA of the fiber, then your analysis should be ok.
Suppose we define the following $\zeta = \ln{\varepsilon}$ and $\xi = \ln{\mu}$, where $\varepsilon$ and $\mu$ are the permittivity and permeability, respectively. In a system with no sources (i.e., $\mathbf{j} = 0$ and $\rho_{c} = 0$), then we know that $\nabla \cdot \mathbf{D} = 0$, where $\mathbf{D} = \varepsilon \ \mathbf{E}$ and $\mathbf{B} = \mu \ \mathbf{H}$. After a little vector calculus we can show that:
$$
\nabla \cdot \mathbf{E} = - \nabla \zeta \cdot \mathbf{E} \tag{0}
$$
Using this and some manipulation of Faraday's law and Ampêre's law, we can show that the general differential equation in terms of electric fields only is given by:
$$
\left( \mu \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \nabla^{2} \right) \mathbf{E} = \left( \mathbf{E} \cdot \nabla \right) \nabla \zeta + \left( \nabla \zeta \cdot \nabla \right) \mathbf{E} + \nabla \left( \zeta + \xi \right) \times \left( \nabla \times \mathbf{E} \right) \tag{1}
$$
We can get a tiny amount of reprieve from this by assuming that the permeability is that of free space, i.e., $\nabla \xi = 0$. If we further argue that the only direction in which gradients matter is along $\hat{x}$ and that the incident wave vector, $\mathbf{k}$, is parallel to this, then we can further simplify Equation 1 to:
$$
\left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) \mathbf{E} = \left( \mathbf{E} \cdot \nabla \right) \zeta' \hat{x} + \left( \zeta' \frac{ \partial }{ \partial x } \right) \mathbf{E} + \zeta' \hat{x} \times \left( \nabla \times \mathbf{E} \right) \tag{2}
$$
where $\zeta' = \tfrac{ \partial \zeta }{ \partial x }$.
After some more manipulation, we can break this up into components to show that:
$$
\begin{align}
\text{x : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{x} & = E_{x} \zeta'' + \zeta' \frac{ \partial E_{x} }{ \partial x } \tag{2a} \\
\text{y : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{y} & = 0 \tag{2b} \\
\text{z : } \left( \mu_{o} \varepsilon \frac{ \partial^{2} }{ \partial t^{2} } - \frac{ \partial^{2} }{ \partial x^{2} } \right) E_{z} & = 0 \tag{2c}
\end{align}
$$
In the limit where the incident wave is entirely transverse, then $E_{x} = 0$ and the x-component (Equation 2a) is entirely zero.
Next you assume that $\mathbf{E} = \mathbf{E}_{o}\left( x \right) e^{i \omega t}$, where $\omega$ is the frequency of the incident wave. Then there will be incident, reflected, and transmitted contributions to the total field at any given point (well the transmitted is always zero in the first medium, of course). Any incident and transmitted contributions with have $\mathbf{k} \cdot \mathbf{x} > 0$ while reflected waves will satisfy $\mathbf{k} \cdot \mathbf{x} < 0$. You define the ratio of the reflected-to-incident fields (well impedances would be more appropriate) to get the coefficient of reflection.
Simpler Approach
A much simpler approach is to know where to look for the answer to these types of questions. As I mentioned in the comments, there has been a ton of work on this very topic (i.e., spatially dependent index of refraction) done for the ionosphere. If we look at, for instance, Roettger [1980] we find a nice, convenient equation for the reflection coefficient, $R$, as a function of the index of refraction, given by:
$$
R = \int \ dx \ \frac{ 1 }{ 2 \ n\left( x \right) } \frac{ \partial n\left( x \right) }{ \partial x } \ e^{-i \ k \ x} \tag{3}
$$
There is no analytical expression for $R$ for your specific index of refraction. However, numerical integration is not difficult if one knows the values for $d$ and $\epsilon$. Note that if we do a Taylor expansion for small $\epsilon$, then the integrand (not including the exponential) is proportional to cosine, to first order in $\epsilon$ (cosine times sine if we go to second order).
References
- Gossard, E.E. "Refractive index variance and its height distribution in different air masses," Radio Sci. 12(1), pp. 89-105, doi:10.1029/RS012i001p00089, 1977.
- Roettger, J. "Reflection and scattering of VHF radar signals from atmospheric refractivity structures," Radio Sci. 15(2), pp. 259-276, doi:10.1029/RS015i002p00259, 1980.
- Roettger, J. and C.H. Liu "Partial reflection and scattering of VHF radar signals from the clear atmosphere," Geophys. Res. Lett. 5(5), pp. 357-360, doi:10.1029/GL005i005p00357, 1978.
Best Answer
The classic, definitive reference for this stuff is:
Snyder & Love, "Optical Waveguide Theory"
My edition is the Chapman and Hall one.
The ratio you seek naturally depends on the fibre's refractive index profile. Snyder and Love give the formula for this ratio for many, many different situations in the book. The two most useful are:
Weak Guidance Step Index Profile Approximation:
The weak guidance approximation to the step-refractive-index profile fibre (S&L, §14.6, p313 in the book I have linked), whose core power fraction is given by:
$$\eta = \frac{W^2}{V^2} + \left(\frac{U\,K_0(W)}{V\,K_1(W)}\right)^2\tag{1}$$
where:
$$U^2+W^2 = V^2\tag{2}$$
and the eigenvalue equation for this fibre is:
$$U\frac{J_1(U)}{J_0(U)}=W\frac{K_1(W)}{K_0(W)}\tag{3}$$
and, as you know:
$$V = k\,\rho\,\sqrt{n_{co}^2-n_{cl}^2}\tag{4}$$
where $\rho$ is core radius, $k$ the freespace wavenumber and $n_{co}, \, n_{cl}$ the core and cladding indices respectively. $J_n$ is the $n^{th}$ order Bessel function of the first kind (the nonsingular one) and $K_n$ the the $n^{th}$ modified Bessel function of the second kind.
What you must do is, given $V$, solve (3) subject to (2) (there will be more than one solution if $V>\omega_{0,1}\approx 2.405$ where $\omega_{0,1}\approx 2.405$ is the first zero $J_0(z)$.
This is all a bit unwieldy and you'll need to build something like a Mathematica notebook to calculate with. However, there is an approximation to this one that works well:
Gaussian Approximation to the Weak Guidance Step Index Profile Approximation
For many if not most practical fibres, the modes look very much like generalised Gaussian modes (in particular Hypergeometric Gaussian modes) whose spotsizes do not vary with propagation (i.e. the optical fibre yields the right amount of focussing to exactly cancel diffraction). One of the major techniques in S&L is the use of variational principles to find best fit Gaussian modes for different fibres. For the step refractive index fibre, the stationary-Lagrangian Gaussian spotsize for the lowest order fibre mode is defined by:
$$\sigma = \frac{\rho}{\sqrt{2\,\log V}}\tag{5}$$
so that the amplitude profile is:
$$E(r) = \exp\left(-\frac{r^2}{2\,\sigma^2}\right) = \exp\left(-\log(V)\,\frac{r^2}{\rho^2}\right)\tag{6}$$
and intensity profile is:
$$I(r) = \exp\left(-2\,\log(V)\,\frac{r^2}{\rho^2}\right)\tag{7}$$
where $\rho$ is the core radius. From these formulas, the fraction of power within a core of radius $\rho$ is:
$$\eta = 1 - \exp\left(-2\,\log(V)\right) = 1-V^{-2}\tag{8}$$
You need to be a little careful with radius and mode field definitions: there are many around. Most often, the quantity called the "mode field diameter" is the Petermann II diameter. For a Gaussian beam, this is the width of the beam wherein the intensity is greater than $1/e^2$ times its peak value, although the Petermann II concept is rather more general than this idea and it has to do with diffracted beamwidths. From (6), the intensity has dropped to $1/e^2$ times the peak value at a radius of $\sqrt{2}\,\sigma$ from the core, so the Petermann II mode field diameter is $2\,\sqrt{2}\,\sigma$. The intensity profile is $I(r) \propto \exp\left(-\frac{r^2}{\sigma^2}\right)$ so that the fraction of power in the core or radius $\rho$ is:
$$\eta = \frac{\int_0^\rho\,r\,\exp\left(-\frac{r^2}{\sigma^2}\right)\,{\rm d}\,r}{\int_0^\infty\,r\,\exp\left(-\frac{r^2}{\sigma^2}\right)\,{\rm d}\,r} = 1-\exp\left(-\frac{\rho^2}{\sigma^2}\right) = 1-\exp\left(-\frac{8\, \rho^2}{\mathscr{P}^2}\right)\tag{9}$$
where $\mathscr{P}$ is the Petermann II diameter or "mode field diameter". So your formula is correct if $w$ stands for the mode field radius NOT diameter. For English speakers, $w$, being the first letter of "width", bespeaks a diameter, so this is confusing. It has to do with the standard notation for the waist parameter of a Gaussian beam, which is the mode field diameter divided by 2. Be careful. the $1-V^{-2}$ formula for the power fraction is much less confusing.
Another formula, very near in its calculated values to (5) for $1.2<V<\infty$ is Marcuse's Petermann II Radius formula:
$$\frac{w}{\rho}=\frac{\mathscr{P}}{2\,\rho}=0.65+\frac{1.619}{V^{\frac{3}{2}}}+\frac{2.879}{V^6}\tag{10}$$
and sometimes this is refined further to:
$$\frac{w}{\rho}=\frac{\mathscr{P}}{2\,\rho}=0.65+\frac{1.619}{V^{\frac{3}{2}}}+\frac{2.879}{V^6}-\left(0.016+\frac{1.561}{V^7}\right)\tag{10}$$
Further Words of Warning
Although the above formulas work for many situations aside from the step index profile if one sticks to $V$ parameters instead of using the core radius directly. That is, the use of the above formulas gives an "effective" core radius for non step index profiles. Mostly, the formulas are extremely insensitive to the exact profiles, so your formulas are adequate. There is one situation in particular where all this can fail spectacularly and that is for dispersion shifted communications fibres. Here the refractive index profile is tailored to shift the zero dispersion point (i.e. wavelength where the derivative of group velocity with respect to wavelength) is in the middle of the 1550nm communications band, whereas it naturally tends to fall around 1300nm for silica fibres. A complicated refractive index profile with several concentric regions of different index is used to achieve the dispersion shift. In this case, you need to check the literature for more appropriate formulas.