Does this generalized second order eigenproblem have exact solution?

I want to study a simple case of Schrodinger equation with position dependent mass and have come up with this. Let:

$$\frac{1}{m(x)}=1+\mu \sin \pi x, \qquad |\mu|<1$$

Now consider a particle with position dependent mass stuck in a one-dimensional box. Solve:

$$-\frac{1}{2} \frac{d}{dx} \left(\frac{1}{m(x)} \frac{d \psi(x)}{dx} \right)=E\psi(x) \\ \psi(0)=\psi(1)=0$$

There's a variety of numerical and approximate methods to solve this, but I wanted to see if there's exact solution for eigenvalues $E_n$ and eigenfunctions $\psi_n$ first, since that would be the best way to check the accuracy of various numerical methods.

Expanding, we obtain:

$$(1+\mu \sin \pi x)~\psi''+\pi \mu \cos \pi x ~\psi' +2E\psi=0$$

Mathematica doesn't give me the general solution, but I'm pretty sure it exists, since I deliberately chose the trigonometric dependence for the reciprocal mass.

My ODEs are quite rusty, so I would appreciate any hint or a reference which would allow me to obtain an exact solution (which may involve special functions and transcendental equation for the eigenvalues of course).

The only ideas I have so far is to try Fourier expansion (which leads to an infinite matrix, so I'm pretty sure would only work as an approximate method), or Fourier transform (involves convolutions, so I'd prefer a more simple method).

If you'd like to know my motivation, I really just want a simple test case for various finite difference schemes, because I'm not sure what's the best way to take into account the position dependence of the mass.

By the way, this is also a good problem for perturbation treatment, since $\mu$ makes for a nice perturbation parameter.

Maple does give a general solution for your differential equation, involving a Heun function

$$y \left( x \right) ={ \left( c_{{1}}{\it HeunG} \left( {\frac {\mu-1}{ 2\,\mu}},{\frac { \left( 5\,\mu-1 \right) {\pi}^{2}+8\,E}{8\,{\pi}^{2} \mu}},{\frac{1}{2}},{\frac{3}{2}},{\frac{3}{2}},{\frac{1}{2}},{\frac { \sin \left( \pi\,x \right) }{2}}+{\frac{1}{2}} \right) \sqrt {-2\,\mu \, \left( \cos \left( \pi\,x \right) \right) ^{2}+2\, \left( \sin \left( \pi\,x \right) +1 \right) \left( 1+\mu \right) }+c_{{2}}{\it HeunG} \left( {\frac {\mu-1}{2\,\mu}},{\frac {E}{{\pi}^{2}\mu}},0,1,{ \frac{1}{2}},{\frac{1}{2}},{\frac {\sin \left( \pi\,x \right) }{2}}+{ \frac{1}{2}} \right) \sqrt {1+\mu\,\sin \left( \pi\,x \right) } \right) {\frac {1}{\sqrt [4]{ \left( \cos \left( \pi\,x \right) \right) ^{2}{\mu}^{2}-2\,\mu\,\sin \left( \pi\,x \right) -{\mu}^{2}-1 }}}} $$

However, I don't think you'll get a nice eigenvalue problem here: note that the solution depends only on $\sin(\pi x)$ and $\cos(\pi x)^2$, so if $c_1$ and $c_2$ are chosen to make it $0$ at $x=0$ then it is also $0$ at $x=1$.

The substitution $u = \sin(\pi x) $ followed by another substitution $v = (1-u)/2 $ leads to a following ODE:

\begin{eqnarray} \psi^{''}_v + \left( \frac{1/2}{v} + \frac{-1/2}{v-1} + \frac{1}{v - \frac{1+\mu}{2 \mu}} \right) \cdot \psi^{'}_v + \frac{-E/(\mu \pi^2)}{v(v-1)(v- \frac{1+\mu}{2 \mu})} \cdot \psi_v = 0 \end{eqnarray}

This is a Heun equation with $a,q = (1+\mu)/(2 \mu),E/(\mu \pi^2)$ and $\alpha,\beta,\gamma,\delta,\epsilon = 0,0,1/2,-1/2,1$.

According to here that equation has two linearly independent solutions given as:

\begin{eqnarray} \psi_v = c_1 \cdot Hl(a,q;\alpha,\beta,\gamma,\delta,v) + c_2 \cdot Hl(a,(a \delta+\epsilon)(1-\gamma)+q,\alpha+1-\gamma,\beta+1-\gamma,2-\gamma,\delta,v) \end{eqnarray}

Clear[phi]; x =.; mu =.; EE =.;
eqn = (1 + mu Sin[Pi x]) D[phi[x], {x, 2}] + 
   Pi mu Cos[Pi x] D[phi[x], {x, 1}] + 2 EE phi[x];

eqn = eqn /. 
      Derivative[1][phi][x] :>  
       Pi Sqrt[1 - u^2] Derivative[1][phi][u] /. 
     Derivative[2][phi][x] :>  
      Pi Sqrt[1 - u^2] D[Pi Sqrt[1 - u^2] Derivative[1][phi][u], 
        u] /. Sin[Pi x] :> u /. Cos[Pi x] :> Sqrt[1 - u^2] /. 
  phi[x] :> phi[u]

eqn = Collect[eqn, Derivative[n_][phi][u], Simplify]

eqn = Collect[
   eqn /. u :> 1 - 2 v /. 
     Derivative[n_.][phi][1 - 2 v] :> Derivative[n][phi][v]/(-2)^n /. 
    phi[1 - 2 v] :> phi[v], Derivative[n_][phi][v], Simplify];

coeff = Coefficient[eqn, Derivative[2][phi][v]];

eqn = Collect[eqn/coeff, Derivative[n_][phi][v], Simplify]
Apart[Coefficient[eqn, Derivative[1][phi][v]]]