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,