checkodesol doesn't recognise solution with unevaluatable integral

[oscarbenjamin]

checkodesol should report True/False to say whether an expression is a solution of an ODE e.g.:

In [89]: eq = f(x).diff(x) - f(x)                                                                                                              

In [90]: eq                                                                                                                                    
Out[90]: 
        d       
-f(x) + ──(f(x))
        dx      

In [91]: sol = exp(-x)                                                                                                                         

In [92]: sol                                                                                                                                   
Out[92]: 
 -x
ℯ  

In [93]: checkodesol(eq, sol)                                                                                                                  
Out[93]: (False, -2)

In [94]: sol = exp(x)                                                                                                                          

In [95]: checkodesol(eq, sol)                                                                                                                  
Out[95]: (True, 0)

However with the following ODE checkodesol gets it wrong:

In [83]: neq = -C1/x - 2*x*f(x) - f(x) + Derivative(f(x), x)                                                                                   

In [84]: neq                                                                                                                                   
Out[84]: 
  C₁                     d       
- ── - 2⋅x⋅f(x) - f(x) + ──(f(x))
  x                      dx      

In [85]: wsol = C2*exp(x**2 + x) + exp(x**2 + x)*Integral(C1*exp(-x**2 - x)/x, x)                                                              

In [86]: wsol                                                                                                                                  
Out[86]: 
                     ⌠                
                     ⎮        2       
     2         2     ⎮     - x  - x   
    x  + x    x  + x ⎮ C₁⋅ℯ           
C₂⋅ℯ       + ℯ      ⋅⎮ ──────────── dx
                     ⎮      x         
                     ⌡                

In [87]: checkodesol(neq, wsol)                                                                                                                
Out[87]: 
(False,
 (2*C1*x*Integral(exp(-x*(x + 1))/x, x) + C1*Integral(exp(-x*(x + 1))/x, x) + 2*C2*x + C2 - (C2 + Integral(C1*exp(-x*(x + 1))/x, x))*(2*x + 1))*exp(x*(x + 1)))

We can correctly find that this is a solution much quicker when using

In [88]: simplify(neq.subs(f(x).diff(x), wsol.diff(x)).subs(f(x), wsol))                                                                       
Out[88]: 0

[jmig5776]

I would like to work on this issue.

[jmig5776]

@oscarbenjamin I had solved it in #15957 Please review.