The neutrino oscillation probability in matter is given as:
where
Now what I do not understand is that "As the energy increases, the probability of oscillation within the sun through the matter effect increases, so the survival probability decreases". I have read this (page 28) in couple of books but unable to cross check it through formula because as I try to increase the energy, both the $sin^2$ decreases.
So what's going wrong here?
I think the crucial point here is that in the sun you have two variables, the neutrino energy $E$ and the electron number density $N_e$. These two enter the matter potential. In order to run into the MSW resonance (and therefore effectively oscillate into muon neutrinos) these two variables have to combine in a way that $\frac{\Delta V}{\Delta m^2} - \cos{2\theta}\approx 0$.
To really check this claim, you would have plug in the best fit values for vacuum mixing for $\Delta m^2$ and $\theta$, but most importantly the continuum of electron number densities of the sun $N_e(x)$.
The transition probability $P_{\nu_e \rightarrow \nu_{\mu}}$ is indeed decreasing with neutrino production energy, i.e., the survival probability $P_{\nu_e \rightarrow \nu_{e}}=1-P_{\nu_e \rightarrow \nu_{\mu}}$ is increasing. Why? The first sine term ($\sin^2 2 \theta_m $) in the oscillation equation (your first equation) determines the amplitude of the oscillations. The limit of that term as $V \rightarrow \infty$, i.e., as energy or density increases towards infinity, is zero: $$ \begin{align} \sin 2 \theta_m = \frac{\sin 2 \theta}{ \sqrt{ \left(\Delta V / \Delta m^2 - \cos 2 \theta \right)^2 + \sin^2 2 \theta} } \\ \lim_{V\to \infty} \; \left[ \frac{\sin 2 \theta}{ \sqrt{ \left( V / \Delta m^2 - \cos 2 \theta \right)^2 + \sin^2 2 \theta} } \right] = 0 \end{align} $$
As $\sin^2 2 \theta_m$ descends towards zero asymptotically with increasing neutrino production energy (and/or production electron number density) above MSW resonance any oscillations in the second sine term, $\sin^2 \left( 1.27 \Delta m_m^2 \frac{L}{E} \right)$, of the oscillation equation will accordingly be reduced in amplitude, going to zero in the limit as well.
The second sine term is not monotonically decreasing (it might appear to be so if only a few samples are calculated), but is rather oscillating rapidly, hence the smeared oscillations in the graph.
The equations you gave should not be used for solar neutrino analysis except at production energies much lower than the MSW resonance.
You can obtain proper equations in the PDG 2018 Review of Particle physics section 14, a free download at lbl.gov.
If the above terse statement of the problem is insufficient you can view an article demonstrating how to graph using the posted equations (with Python code) and how they malfunction at higher energies (with a link to a proper equation and an article describing that) at this link: MSW flavor calculation
