What's wrong with :
 sage: E, k=var("E, k")
sage: y=function("y")
sage: de=y(x).diff(x, 2)*x+y(x).diff(x)*(2*k+1-2*x^2)+y(x)*(E-2*k)*x==0
sage: desolve(de,y(x),ivar=x, contrib_ode=True)
[y(x) == _K1*hypergeometric_M(-1/4*E + 1/2*k, k + 1, x^2) + _K2*hypergeometric_U(-1/4*E + 1/2*k, k + 1, x^2)]
 [ FWIW, Sympy's dsolve doesn't seem to be able to solve this one...]
 Neither Sage, Fricas, Giac nor Sympy seem to be able to check "easily" this soluion analytically. But Mathematica (the gratis WolframEngine 13.2 here) can :
 sage: Sol=desolve(de,y(x),ivar=x, contrib_ode=True) ; Sol # Same as before...
[y(x) == _K1*hypergeometric_M(-1/4*E + 1/2*k, k + 1, x^2) + _K2*hypergeometric_U(-1/4*E + 1/2*k, k + 1, x^2)]
sage: y0(x)=Sol[0].rhs()
sage: de.substitute_function(y,y0)._mathematica_().FullSimplify()
True
 Alternatively :
 sage: de.substitute_function(y,y0).lhs()._mathematica_().FullSimplify()
0
 A numerical check could be done with Sage alone by substituting random numerical values for _K1 and _K2 and plotting the de's left hand on the domain(s) of interest.
 HTH,