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$$.

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.