Monotonicity of Schrödinger Eigenvalues Let us consider the Schrödinger operator
$$
H_hf(x)=-\frac{d^2}{dx^2}f(x)+h(h\sin^2(x)-\cos(x))f(x)
$$
on $L^2[-\pi,\pi]$ with Neumann boundary conditions $f^\prime(\pm\pi)=0$. Here, $h\geq 0$ is a parameter. 
It is easy to see that (up to normalization) 
$$
\psi_0^h(x)=e^{h\cos(x)}
$$
is an eigenfunction for the eigenvalue $0$. 
Running some numerics suggests that the mapping $h\mapsto\lambda_n(H_h)$ where $\lambda_n$ denotes the $n^\text{th}$ eigenvalue of $H_h$, $n=0,1,2,\dots$, is monotone non-decreasing. The question is how to prove this.
My approach was the following: By a famous result,
$$
\frac{d}{dh}\lambda_n(H_h)=\langle \psi_n^h,H_h^\prime\psi_n^h\rangle
$$
where $\psi_n^h$ is the $n^\text{th}$ eigenfuntion of $H_h$ and $H_h^\prime f(x)=(2h\sin^2(x)-\cos(x))f(x)$. However, I am struggling to prove that $\langle \psi_n^h,H_h^\prime\psi_n^h\rangle\geq 0$.
Update
Maybe we should first focus on the second eigenvalue $\lambda_1(H_h)$. By general theory, $\psi_1^h$ is a differentiable odd function. Furthermore, $\psi_1^h$ has only one zero (at $x=0$). To obtain $\frac{d}{dh}\lambda_n(H_h)\geq 0$, it would thus be sufficient to prove
$$
\max_{x\in[0,x_0]}|\psi_1^h(x)|\leq\min_{x\in[x_0,\pi]}|\psi_1^h(x)|,
$$
where $x_0=2\arctan\left(\sqrt{\sqrt{16h^2 + 1} - 4 h}\right)$ is the zero of $2h\sin^2(x)-\cos(x)$ in $[0,\pi]$, since
$$
\int_{0}^\pi 2h\sin^2(x)-\cos(x)\,dx=h\pi\geq 0.
$$
Update 2
This is based on the comment by Carlo Beenakker. We transform the initial Neumann problem into a Dirichlet problem. A computation gives that the non-zero spectrum of $H_h$ coincides with the spectrum of
\begin{align}
K_hf(x)&=-f''(x)+h(h\sin^2(x)+\cos(x))f(x)\\&=-f''(x)+\left(h\cos(x)-\frac{h^2}{2}\cos(2x)+\frac{h^2}{2}\right)f(x)
\end{align} 
acting on $L^2[-\pi,\pi]$ with Dirichlet boundary conditions $f(\pm\pi)=0$.
 A: As requested in the comment, let me explain what I could extract from the literature on this problem. (This is not the solution asked for in the OP, but what I have would be too long for a comment.)
The potential in the OP is known as the "Razavy potential" or "double cosine potential". A recent study is Exact Solutions of the Razavy Cosine Type Potential (2018). The Schrödinger equation is
$$-\psi''(x)+V(x)\psi(x)=E\psi(x),\;\;V(x)=\tfrac{1}{4}\xi^2\sin^2 x-(a+1)\xi\cos x,$$
on $-\pi\leq x\leq\pi$ with $\psi(\pm\pi)=0$. The differential equation in the OP is for $\xi=2h$ and $a=-3/2$. (I note that the cited paper assumes $\xi,a>0$, but that does not seem to be an essential condition on the solution.)
The substitution $\psi(x)=\exp(\tfrac{1}{2}\xi\cos x)\phi(x)$ and the change of variables $z=\cos^2(x/2)$ produces a confluent Heun differential equation with solution given by the Heun function:
$$\psi(x)=\exp(\tfrac{1}{2}\xi\cos x)H\bigl(2\xi,-1/2,-1/2,-(2a+1)\xi,2(a+1)\xi+3/8-E;\cos^2(x/2)\bigr).$$
The energy $E$ should then be obtained from the boundary condition $x=\pm\pi$, so at the origin for the Heun function, but the cited paper does not succeed in obtaining a closed-form solution.
