Repository: sympy/sympy
Pull request: Fixed checkodesol function
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234". See
https://github.com/blog/1506-closing-issues-via-pull-requests . Please also
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #15913 

#### Brief description of what is fixed or changed
Earlier
```
eq = -C1/x - 2*x*f(x) - f(x) + Derivative(f(x), x)
sol = C2*exp(x**2 + x) + exp(x**2 + x)*Integral(C1*exp(-x**2 - x)/x, x)
checkodesol(eq, sol)
(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)))
```
Now it returns correct answer
```
eq = -C1/x - 2*x*f(x) - f(x) + Derivative(f(x), x)
sol = C2*exp(x**2 + x) + exp(x**2 + x)*Integral(C1*exp(-x**2 - x)/x, x)
checkodesol(eq, sol)
(True, 0)
```

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* solvers
  * checkodesol is more robust when checking solutions involving Integrals
<!-- END RELEASE NOTES -->


Problem statement:
checkodesol doesn't recognise solution with unevaluatable integral
checkodesol should report True/False to say whether an expression is a solution of an ODE e.g.:
```julia
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:
```julia
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
```julia
In [88]: simplify(neq.subs(f(x).diff(x), wsol.diff(x)).subs(f(x), wsol))                                                                       
Out[88]: 0
```

Hints from issue discussion:
I would like to work on this issue.
@oscarbenjamin I had solved it in #15957 Please review.

Commit message: Update ode.py

Patch to review:
diff --git a/sympy/solvers/ode.py b/sympy/solvers/ode.py
index 8fc6f514a02b..bae3cb0eb974 100644
--- a/sympy/solvers/ode.py
+++ b/sympy/solvers/ode.py
@@ -338,14 +338,13 @@
     )
 
 
-def sub_func_doit(eq, func, new):
+def sub_func_doit(eq, func, new, order=None):
     r"""
     When replacing the func with something else, we usually want the
     derivative evaluated, so this function helps in making that happen.
-
-    To keep subs from having to look through all derivatives, we mask them off
-    with dummy variables, do the func sub, and then replace masked-off
-    derivatives with their doit values.
+    If ``order`` is provided, then all derivatives of all orders will
+    be evaluated; otherwise, only the derivatives present (default)
+    will be evaluated.
 
     Examples
     ========
@@ -363,25 +362,15 @@ def sub_func_doit(eq, func, new):
     x*(-1/(x**2*(z + 1/x)) + 1/(x**3*(z + 1/x)**2)) + 1/(x*(z + 1/x))
     ...- 1/(x**2*(z + 1/x)**2)
     """
-    reps = {}
-    repu = {}
-    for d in eq.atoms(Derivative):
-        u = Dummy('u')
-        repu[u] = d.subs(func, new).doit()
-        reps[d] = u
-
-    # Make sure that expressions such as ``Derivative(f(x), (x, 2))`` get
-    # replaced before ``Derivative(f(x), x)``:
-    #
-    # Also replace e.g. Derivative(x*Derivative(f(x), x), x) before
-    # Derivative(f(x), x)
-    def cmp(subs1, subs2):
-        return subs2[0].has(subs1[0]) - subs1[0].has(subs2[0])
-    key = lambda x: (-x[0].derivative_count, cmp_to_key(cmp)(x))
-    reps = sorted(reps.items(), key=key)
-
-    return eq.subs(reps).subs(func, new).subs(repu)
-
+    reps = {func: new}
+    if order:
+        x = func.args[0]
+        for i in range(order, -1, -1):
+            reps[func.diff(x, i)] = new.diff(x, i)
+    else:
+        for d in eq.atoms(Derivative):
+            reps[d] = d.subs(func, new).doit(deep=False)
+    return eq.xreplace(reps)
 
 def get_numbered_constants(eq, num=1, start=1, prefix='C'):
     """
@@ -2415,7 +2404,7 @@ def checkodesol(ode, sol, func=None, order='auto', solve_for_func=True):
             ode_diff = ode.lhs - ode.rhs
 
             if sol.lhs == func:
-                s = sub_func_doit(ode_diff, func, sol.rhs)
+                s = sub_func_doit(ode_diff, func, sol.rhs, order)
             else:
                 testnum += 1
                 continue
diff --git a/sympy/solvers/tests/test_ode.py b/sympy/solvers/tests/test_ode.py
index 59d9794dec77..acfe25ed7f35 100644
--- a/sympy/solvers/tests/test_ode.py
+++ b/sympy/solvers/tests/test_ode.py
@@ -3115,3 +3115,8 @@ def test_factoring_ode():
     soln = Eq(f(x), (C1*x**2/2 + C2*x + C3 - x)/(1 + x))
     assert checkodesol(eqn, soln, order=2, solve_for_func=False)[0]
     assert soln == dsolve(eqn, f(x))
+
+def test_issue_15913():
+    eq = -C1/x - 2*x*f(x) - f(x) + Derivative(f(x), x)
+    sol = C2*exp(x**2 + x) + exp(x**2 + x)*Integral(C1*exp(-x**2 - x)/x, x)
+    assert checkodesol(eq, sol) == (True, 0)


Code review comments:
diff_hunk: @@ -363,25 +362,15 @@ def sub_func_doit(eq, func, new):
     x*(-1/(x**2*(z + 1/x)) + 1/(x**3*(z + 1/x)**2)) + 1/(x*(z + 1/x))
     ...- 1/(x**2*(z + 1/x)**2)
     """
-    reps = {}
-    repu = {}
-    for d in eq.atoms(Derivative):
-        u = Dummy('u')
-        repu[u] = d.subs(func, new).doit()
-        reps[d] = u
-
-    # Make sure that expressions such as ``Derivative(f(x), (x, 2))`` get
-    # replaced before ``Derivative(f(x), x)``:
-    #
-    # Also replace e.g. Derivative(x*Derivative(f(x), x), x) before
-    # Derivative(f(x), x)
-    def cmp(subs1, subs2):
-        return subs2[0].has(subs1[0]) - subs1[0].has(subs2[0])
-    key = lambda x: (-x[0].derivative_count, cmp_to_key(cmp)(x))
-    reps = sorted(reps.items(), key=key)
-
-    return eq.subs(reps).subs(func, new).subs(repu)
-
+    reps = {func: new}
+    if order:
+        x = func.args[0]
+        for i in range(order, -1, -1):
+            reps[func.diff(x, i)] = new.diff(x, i)
+    else:
+        for d in eq.atoms(Derivative):
+            reps[d] = d.subs(func, new).doit(deep=False)
path: sympy/solvers/ode.py
text: @user1:
the only other thing that might make this slow in comparison to the above is that this is using `d.subs` rather than computing the derivative directly...maybe the subs is slow and if the rhs. But now another complicating factor occurs to me: any given derivative might not be constant-free, e.g. `Derivative(y*f(x),x).subs(f(x).diff(x),x)` in which case the above would fail.

So perhaps:
```python
reps= {func: new}
for d in eq.atoms(Derivative):
    if d.expr == func:
        reps[d] = new.diff(*d.variable_count)
    else:
        reps[d] = d.xreplace({func: new}).doit(deep=False)
return eq.xreplace(reps)
```

@user1:
a test of that complicating issue should be added

@user1:
This is not necessary since `_preprocess` already makes sure derivatives are evaluated.