Repository: sympy/sympy
Pull request: [WIP] Rational Riccati solver
This PR introduces a Rational Riccati ODE solver. This code can also be used to solve a subclass of 2nd order ODEs solvable by [Kovacic's Algorithm](https://cs.uwaterloo.ca/research/tr/1984/CS-84-35.pdf)

#### References to other Issues or PRs
Partial Fix #19183, #21503

#### Brief description of what is fixed or changed
The algorithm present [here](https://www3.risc.jku.at/publications/download/risc_5387/PhDThesisThieu.pdf) is implemented. The examples used to test are [here](https://www3.risc.jku.at/publications/download/risc_5197/RISCReport15-19.pdf).

#### Other comments
I will soon add the Kovacic solver.

#### Release Notes
<!-- BEGIN RELEASE NOTES -->
* solvers
  * ode:  Added a new solver for Riccati equation
<!-- END RELEASE NOTES -->

Problem statement:
Riccati's equation
Hello
I have a question regarding using Sympy for solving Riccati's equation:
y’  = -y**2  + x*y + 1    
It seems that Sympy is not able to solve this equation, while a closed-formed solution in terms of the error function exists (Mathematica has no problem finding this solution).
Would someone have any advice?
Thanks a lot for your time.
Stay safe,
Elie.

Hints from issue discussion:
Sympy currently gives a series solution:
```
In [3]: dsolve(Eq(f(x).diff(x), -f(x)**2 + f(x)*x + 1))                                                                           
Out[3]: 
                      3 ⎛      2   ⎛      2⎞ ⎛    2    ⎞    ⎞    5 ⎛    2 ⎛        2⎞        2   ⎛      2⎞ ⎛     2    ⎞   ⎛     
         ⎛      2⎞   x ⋅⎝- 3⋅C₁  + ⎝1 - C₁ ⎠⋅⎝6⋅C₁  - 1⎠ + 1⎠   x ⋅⎝8⋅C₁ ⋅⎝1 - 3⋅C₁ ⎠ - 14⋅C₁  + ⎝1 - C₁ ⎠⋅⎝36⋅C₁  - 5⎠ + ⎝1 - C
f(x) = x⋅⎝1 - C₁ ⎠ + ──────────────────────────────────────── + ────────────────────────────────────────────────────────────────
                                        6                                                                                       

 2⎞ ⎛       2 ⎛        2⎞        2     ⎛        2⎞ ⎛      2⎞    ⎞    ⎞            2 ⎛    2    ⎞       4 ⎛     2     ⎛        2⎞ 
₁ ⎠⋅⎝- 16⋅C₁ ⋅⎝1 - 3⋅C₁ ⎠ + 36⋅C₁  + 8⋅⎝1 - 9⋅C₁ ⎠⋅⎝1 - C₁ ⎠ - 5⎠ + 2⎠        C₁⋅x ⋅⎝2⋅C₁  - 1⎠   C₁⋅x ⋅⎝12⋅C₁  + 8⋅⎝1 - 3⋅C₁ ⎠⋅
────────────────────────────────────────────────────────────────────── + C₁ + ───────────────── + ──────────────────────────────
 120                                                                                  2                                24       

⎛      2⎞    ⎞        
⎝1 - C₁ ⎠ - 5⎠    ⎛ 6⎞
────────────── + O⎝x ⎠
```
There is a solver for one type of Ricatti equation but maybe this is a different one.
Thanks. But a solution in terms of a series is not really useful when a closed-form solution exists
y(x) = E^(-x^2/2)/ (Sqrt(Pi/2) Erf(x/Sqrt(2)) + c)) + x
(with Erf the error function, and c an arbitrary constant).
> a solution in terms of a series is not really useful

I agree. In fact I don't think dsolve should return series solutions at all (#15916).

There is a handler for one type of Riccati ODE. I'm not sure if should handle this case or if another solver should be added.
The solver 'Riccati_special_minus2' is of no help in our case. Maybe there is another one, but I was not able to find it.

This is a special type of Riccati ODE, and it can be solved by finding a rational transformation to a Bernoulli ODE. This has not been implemented in sympy as far as I know. The Lie group is [0,exp((x^2)/2)*(y-x)^2] and this group can also not be found with the current implementations.
Another alternative is transforming it to a second order linear ode and use kovacic algorithm to solve second order linear ODEs with Liouvillian solutions. But it seems that this method is also not implemented in sympy? You could try the jupyter implementation of the maxima cas, the Kovacic algorithm as well as Lie's algorithm work for this ODE.   
Interesting. How exactly would you transform this to a linear 2nd order ODE?

I had a quick look at the Kovacic algorithm and that definitely seems like something that should be added to sympy:
https://www.sciencedirect.com/science/article/pii/S0747717186800104
Q. How exactly would you transform this to a linear 2nd order ODE?

A. Let’s start from the Riccati’s equation (see https://en.wikipedia.org/wiki/Riccati_equation)

 $y'(x) = a(x) + b(x) y(x) + c(x) y^2(x)$ 

with $a(x) \not= 0$ and $c(x) \not= 0$. (Note: If $a(x) = 0$ the equation reduces to a Bernoulli equation, while if $c(x) = 0$ the equation becomes a 1st oder linear ODE)

Let us define $S(x) = a(x) c(x)$ and R(x) = b(x) + c'(x)/c(x)$. 

One can show (see https://en.wikipedia.org/wiki/Riccati_equation) that if $u(x)$ is a solution of the 2nd order linear ODE:

$u''(x) - R(x) u'(x) + S(x) u(x) = 0$,

then $y(x) = -u'(x)/(c(x) u(x)$ is a solution of the original Riccati equation.


> Le 27 avr. 2020 à 02:31, Oscar Benjamin <notifications@github.com> a écrit :
> 
> 
> Interesting. How exactly would you transform this to a linear 2nd order ODE?
> 
> I had a quick look at the Kovacic algorithm and that definitely seems like something that should be added to sympy:
> https://www.sciencedirect.com/science/article/pii/S0747717186800104 <https://www.sciencedirect.com/science/article/pii/S0747717186800104>
> —
> You are receiving this because you authored the thread.
> Reply to this email directly, view it on GitHub <https://github.com/sympy/sympy/issues/19183#issuecomment-619650042>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/APKF3JCQOFGONBZOKAKZWATROTG5ZANCNFSM4MP3R4WQ>.
> 


Thanks so it's basically the substitution `y = u'/cu` i.e.:
```julia
In [19]: a, b, c, y, u = symbols('a, b, c, y, u', cls=Function)                                                                                

In [20]: x = Symbol('x')                                                                                                                       

In [21]: eq_y = Eq(Derivative(y(x), x), a(x) + b(x)*y(x) + c(x)*y(x)**2)                                                                       

In [22]: eq_y                                                                                                                                  
Out[22]: 
d                                   2   
──(y(x)) = a(x) + b(x)⋅y(x) + c(x)⋅y (x)
dx                                      

In [23]: eq_u_subs = eq_y.subs(y(x), -Derivative(u(x), x)/(c(x)*u(x)))                                                                         

In [24]: eq_u_subs                                                                                                                             
Out[24]: 
                                                  2
  ⎛ d        ⎞               d          ⎛d       ⎞ 
  ⎜-──(u(x)) ⎟          b(x)⋅──(u(x))   ⎜──(u(x))⎟ 
d ⎜ dx       ⎟               dx         ⎝dx      ⎠ 
──⎜──────────⎟ = a(x) - ───────────── + ───────────
dx⎝c(x)⋅u(x) ⎠            c(x)⋅u(x)            2   
                                         c(x)⋅u (x)

In [25]: eq_u = Eq(Derivative(u(x), x, 2), solve(eq_u_subs.doit(), Derivative(u(x), x, 2))[0].collect(u(x)))                                   

In [26]: eq_u                                                                                                                                  
Out[26]: 
            ⎛       d       ⎞                          
  2         ⎜       ──(c(x))⎟                          
 d          ⎜       dx      ⎟ d                        
───(u(x)) = ⎜b(x) + ────────⎟⋅──(u(x)) - a(x)⋅c(x)⋅u(x)
  2         ⎝         c(x)  ⎠ dx                       
dx                                                     
```
Yes, your calculation is correct.
(Note that there is a small typo in the first line of your message: the substitution is y = - u’/(c u)  )


> Le 27 avr. 2020 à 11:58, Oscar Benjamin <notifications@github.com> a écrit :
> 
> 
> Thanks so it's basically the substitution y = u'/cu i.e.:
> 
> In [19]: a, b, c, y, u = symbols('a, b, c, y, u', cls=Function)                                                                                
> 
> In [20]: x = Symbol('x')                                                                                                                       
> 
> In [21]: eq_y = Eq(Derivative(y(x), x), a(x) + b(x)*y(x) + c(x)*y(x)**2)                                                                       
> 
> In [22]: eq_y                                                                                                                                  
> Out[22]: 
> d                                   2   
> ──(y(x)) = a(x) + b(x)⋅y(x) + c(x)⋅y (x)
> dx                                      
> 
> In [23]: eq_u_subs = eq_y.subs(y(x), -Derivative(u(x), x)/(c(x)*u(x)))                                                                         
> 
> In [24]: eq_u_subs                                                                                                                             
> Out[24]: 
>                                                   2
>   ⎛ d        ⎞               d          ⎛d       ⎞ 
>   ⎜-──(u(x)) ⎟          b(x)⋅──(u(x))   ⎜──(u(x))⎟ 
> d ⎜ dx       ⎟               dx         ⎝dx      ⎠ 
> ──⎜──────────⎟ = a(x) - ───────────── + ───────────
> dx⎝c(x)⋅u(x) ⎠            c(x)⋅u(x)            2   
>                                          c(x)⋅u (x)
> 
> In [25]: eq_u = Eq(Derivative(u(x), x, 2), solve(eq_u_subs.doit(), Derivative(u(x), x, 2))[0].collect(u(x)))                                   
> 
> In [26]: eq_u                                                                                                                                  
> Out[26]: 
>             ⎛       d       ⎞                          
>   2         ⎜       ──(c(x))⎟                          
>  d          ⎜       dx      ⎟ d                        
> ───(u(x)) = ⎜b(x) + ────────⎟⋅──(u(x)) - a(x)⋅c(x)⋅u(x)
>   2         ⎝         c(x)  ⎠ dx                       
> dx                                                     
> —
> You are receiving this because you authored the thread.
> Reply to this email directly, view it on GitHub <https://github.com/sympy/sympy/issues/19183#issuecomment-619872552>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/APKF3JFSGBWW76462BXP7DDROVJLDANCNFSM4MP3R4WQ>.
>

Commit message: Added riccati submodule

Patch to review:
diff --git a/sympy/solvers/ode/ode.py b/sympy/solvers/ode/ode.py
index 98bad219834d..d97d155371fe 100644
--- a/sympy/solvers/ode/ode.py
+++ b/sympy/solvers/ode/ode.py
@@ -332,6 +332,7 @@
     "Liouville_Integral",
     "2nd_nonlinear_autonomous_conserved",
     "2nd_nonlinear_autonomous_conserved_Integral",
+    "1st_rational_riccati",
     )
 
 lie_heuristics = (
@@ -1053,6 +1054,7 @@ class in it.  Note that a hint may do this anyway if
         Factorable: ('factorable',),
         RiccatiSpecial: ('Riccati_special_minus2',),
         SecondNonlinearAutonomousConserved: ('2nd_nonlinear_autonomous_conserved',),
+        RationalRiccati: ('1st_rational_riccati',),
     }
     for solvercls in solvers:
         solver = solvercls(ode)
@@ -7096,3 +7098,4 @@ def _nonlinear_3eq_order1_type5(x, y, z, t, eq):
 from .single import (NthAlgebraic, Factorable, FirstLinear, AlmostLinear,
         Bernoulli, SingleODEProblem, SingleODESolver, RiccatiSpecial,
         SecondNonlinearAutonomousConserved)
+from .riccati import RationalRiccati
diff --git a/sympy/solvers/ode/riccati.py b/sympy/solvers/ode/riccati.py
new file mode 100644
index 000000000000..dd8cba72e0e4
--- /dev/null
+++ b/sympy/solvers/ode/riccati.py
@@ -0,0 +1,500 @@
+r"""
+This module contains :py:meth:`~sympy.solvers.ode.riccati.RationalRiccati`,
+a solver class which gives rational solutions to first order Riccati ODEs.
+A general first order Riccati ODE is given by
+
+.. math:: y' = b_0(x) + b_1(x)w + b_2(x)w^2
+
+where `b_0, b_1` and `b_2` can be arbitrary functions of `x` with `b_2 \ne 0`.
+
+Background
+==========
+
+A Riccati equation can be transformed to its normal form
+
+.. math:: y' + y^2 = a(x)
+
+using the transformation
+
+.. math:: y = -b_2(x) - \frac{b'_2(x)}{2 b_2(x)} - \frac{b_1(x)}{2}
+
+where `a(x)` is given by
+
+.. math:: a(x) = \frac{1}{4}\left(\frac{b_2'}{b_2} + b_1\right)^2 - \frac{1}{2}\left(\frac{b_2'}{b_2} + b_1\right)' - b_0 b_2
+
+Thus, we can develop an algorithm to solve for the Riccati equation
+in its normal form, which would in turn give us the solution for
+the original Riccati equation.
+
+Algorithm
+=========
+
+The algorithm implemented here was presented in a Ph.D thesis by
+N. Thieu Vo. The entire thesis can be found here -
+https://www3.risc.jku.at/publications/download/risc_5387/PhDThesisThieu.pdf
+
+We have only implemented the Rational Riccati solver (Pg 78 in Thesis).
+Before we proceed towards the implementation of the algorithm, a
+few definitions to understand -
+
+1. Valuation of a Rational Function at `\infty`:
+    The valuation of a rational function `p(x)` at `\infty` is equal
+    to the difference between the degree of the denominator and the
+    numerator of `p(x)`.
+
+    Note: A general definition of valuation of a rational function
+    at any value of `x` can be found in Pg 63 of the thesis, but
+    is not of any interest for this algorithm.
+
+2. Zeros and Poles of a Rational Function:
+    Let `a(x) = \frac{S(x)}{T(x)}, T \ne 0` be a rational function
+    of `x`. Then -
+
+    a. The Zeros of `a(x)` are the roots of `S(x)`.
+    b. The Poles of `a(x)` are the roots of `T(x)`. However, `\infty`
+    can also be a pole of a(x). We say that `a(x)` has a pole at
+    `\infty` if `a(\frac{1}{x})` has a pole at 0.
+
+Every pole is associated with an order that is equal to the multiplicity
+of its appearence as a root of `T(x)`. A pole is called a simple pole if
+it has an order 1. Similarly, a pole is called a multiple pole if it has
+an order `\ge` 2.
+
+Necessary Conditions
+====================
+
+For a Riccati equation in its normal form,
+
+.. math:: y' + y^2 = a(x)
+
+we can define
+
+a. A pole is called a movable pole if it is a pole of `y(x)` and is not
+a pole of `a(x)`.
+b. Similarly, a pole is called a non-movable pole if it is a pole of both
+`y(x)` and `a(x)`.
+
+Then, the algorithm states that a rational solution exists only if -
+
+a. Every pole of `a(x)` must be either a simple pole or a multiple pole
+of even order.
+b. The valuation of `a(x)` at `\infty` must be even or be `\ge` 2.
+
+Implementation
+==============
+
+The code is written to match the steps given in the thesis (Pg 82)
+
+Step 0 : Match the equation -
+Find `b_0, b_1` and `b_2`. If `b_2 = 0` or no such functions exist, raise
+an error
+
+Step 1 : Transform the equation to its normal form as explained in the
+theory section.
+
+Step 2 : Initialize 2 empty set of solutions, ``sol`` and ``presol``.
+
+Step 3 : If `a(x) = 0`, append `\frac{1}/{(x - C1)}` to ``presol``.
+
+Step 4 : If `a(x)` is a rational non-zero number, append `\pm \sqrt{a}`
+to ``presol``.
+
+Step 5 : Find the poles and their multiplicities of `a(x)` using
+``find_poles``. Let the number of poles be `n`. Also find the
+valuation of `a(x)` at `\infty` using ``val_at_inf``.
+
+Step 6 : Find `n` c-vectors (one for each pole) and 1 d-vector using
+``construct_c`` and ``construct_d``. Now, determine all the ``2**(n + 1)``
+combinations of choosing between 2 choices for each of the `n` c-vectors
+and 1 d-vector.
+
+For each of these above combinations, do
+
+Step 8 : Compute `m` in ``compute_degree``
+
+Step 9 : In ``compute_degree``, compute ybar as well. If ybar contains any
+term like `\frac{c_{ij}}{(x - oo)^j}` where `c_{ij}` is non-zero, discard
+this solution.
+
+Step 10 : If `m` is a non-negative integer -
+
+Step 11: Find a polynomial solution of degree `m` for the auxiliary equation.
+
+There are 2 cases possible -
+
+    a. `m` is a non-negative integer: We can solve for the coefficients
+    in `p(x)` using Undetermined Coefficients.
+
+    b. `m` is an expression: In this case, we may not be able to determine
+    if it is possible for `m` to be a non-negative integer. However, the
+    solution to the auxiliary equation always seems to be x**m + C1 (where
+    C1 is an arbitrary constant). This has been confirmed by testing many
+    examples. NOTE that this is NOT mentioned in the thesis and is an
+    assumption made by me. If any failing test cases are found, this case
+    must be removed.
+
+Step 12 : For each `p(x)` that exists, append `ybar + \frac{p'(x)}{p(x)}`
+to ``presol``.
+
+Step 13 : For each solution in ``presol``, apply an inverse transformation
+and append them to ``sol``, so that the solutions of the original equation
+are found using the solutions of the equation in its normal form.
+"""
+
+from .single import SinglePatternODESolver
+
+from sympy.core import S
+from sympy.core.add import Add
+from sympy.core.function import Function
+from sympy.core.numbers import oo
+from sympy.core.relational import Eq
+from sympy.core.symbol import Symbol, Wild
+from sympy.functions import sqrt
+from sympy.polys.polytools import cancel, degree
+from sympy.polys.polyroots import roots
+from sympy.simplify.radsimp import numer, denom
+from sympy.series.limits import limit
+from sympy.solvers.solvers import solve
+
+class RationalRiccati(SinglePatternODESolver):
+    r"""
+    Gives all rational solutions to the first
+    order Riccati Differential Equation
+
+    .. math :: y' = b_0(x) + b_1(x) y + b_2(x) y^2
+
+    where `b_0`, `b_1` and `b_2` are rational functions of `x`
+    with `b_2 \ne 0` (`b_2 = 0` would make it a Bernoulli equation).
+
+    Examples
+    ========
+
+    >>> from sympy import Symbol, Function, dsolve, checkodesol
+    >>> f = Function('f')
+    >>> x = Symbol('x')
+
+    >>> eq = x**3*f(x).diff(x) - x**2*f(x) - f(x)**2
+    >>> sols = dsolve(eq, hint="1st_rational_riccati")
+    >>> sols
+    Eq(f(x), x**2*(-x/(C1 + x) + 1))
+    >>> checkodesol(eq, sols)
+    (True, 0)
+
+    >>> eq = -x**4*f(x)**2 + x**3*f(x).diff(x) + x**2*f(x) + 20
+    >>> sols = dsolve(eq, hint="1st_rational_riccati")
+    >>> sols
+    [Eq(f(x), (-9*x**8/(C1 + x**9) + 4/x)/x), Eq(f(x), 4/x**2)]
+    >>> checkodesol(eq, sols)
+    [(True, 0), (True, 0)]
+
+    >>> eq = x**3*f(x).diff(x) - x**2*(f(x) + f(x).diff(x)) + 2*x*f(x) - f(x)**2
+    >>> sols = dsolve(eq, hint="1st_rational_riccati")
+    >>> sols
+    Eq(f(x), x**2*(1 - (x - 1)/(C1 + x)))
+    >>> checkodesol(eq, sols)
+    (True, 0)
+
+    References
+    ==========
+
+    .. [1] Algorithm (Pg 78) - https://www3.risc.jku.at/publications/download/risc_5387/PhDThesisThieu.pdf
+    .. [2] Examples - https://www3.risc.jku.at/publications/download/risc_5197/RISCReport15-19.pdf
+    """
+    has_integral = False
+    hint = "1st_rational_riccati"
+    order = [1]
+
+    def _wilds(self, f, x, order):
+        b0 = Wild('b0', exclude=[f(x), f(x).diff(x)])
+        b1 = Wild('b1', exclude=[f(x), f(x).diff(x)])
+        b2 = Wild('b2', exclude=[f(x), f(x).diff(x)])
+        return (b0, b1, b2)
+
+    def _equation(self, fx, x, order):
+        b0, b1, b2 = self.wilds()
+        return fx.diff(x) - b0 - b1*fx - b2*fx**2
+
+    def _matches(self):
+        eq = self.ode_problem.eq_expanded
+        f = self.ode_problem.func.func
+        x = self.ode_problem.sym
+        order = self.ode_problem.order
+
+        if order != 1:
+            return False
+
+        eq = eq.collect(f(x))
+        cf = eq.coeff(f(x).diff(x))
+
+        # There must be an f(x).diff(x) term.
+        # eq must be an Add object since we are using the expanded
+        # equation and it must have atleast 2 terms (b2 != 0)
+        if cf != 0 and isinstance(eq, Add):
+
+            # Divide all coefficients by the coefficient of f(x).diff(x)
+            # and add the terms again to get the same equation
+            eq = Add(*map(lambda x: cancel(x/cf), eq.args)).collect(f(x))
+
+            # Get the pattern for the Riccati equation
+            pattern = self._equation(f(x), x, order)
+
+            # Match the equation with the pattern
+            self._wilds_match = match = eq.match(pattern)
+
+            # If there is a match, verify it
+            if match is not None:
+                return self._verify(f(x))
+        return False
+
+    def _verify(self, fx):
+        b0, b1, b2 = self.wilds_match()
+        x = self.ode_problem.sym
+
+        # b0, b1 and b2 must all be rational functions of x, with b2 != 0
+        return all([b2 != 0, b0.is_rational_function(x), b1.is_rational_function(x), b2.is_rational_function(x)])
+
+    def solve_aux_eq(self, a, x, m, ybar):
+        p = Function('p')
+        # Eq (5.16) in Thesis - Pg 81
+        auxeq = numer((p(x).diff(x, 2) + 2*ybar*p(x).diff(x) + (ybar.diff(x) + ybar**2 - a)*p(x)).together())
+
+        # Assume that the solution is of the type
+        # p(x) = C0 + C1*x + ... + Cm*x**m
+        psyms = get_numbered_constants(auxeq, m)
+        if type(psyms) != tuple:
+            psyms = (psyms, )
+        psol = x**m
+        for i in range(m):
+            psol += psyms[i]*x**i
+        if m != 0:
+            # m is a non-zero integer. Find the constant terms using undetermined coefficients
+            return psol, solve(auxeq.subs(p(x), psol).doit().expand(), psyms, dict=True), True
+        else:
+            # m ==0 . Check if 1 (x**0) is a solution to the auxiliary equation
+            cf = auxeq.subs(p(x), psol).doit().expand()
+            return S(1), cf, cf == 0
+
+    def compute_degree(self, x, poles, choice, c, d, N):
+        ybar = 0
+        m = S(d[choice[0]][-1])
+
+        # Calculate the first (nested) summation for ybar
+        # as given in Step 9 of the Thesis (Pg 82)
+        for i in range(len(poles)):
+            for j in range(len(c[i][choice[i + 1]])):
+                # If one of the poles is infinity and its coefficient is
+                # not zero, the given solution is invalid as there will be
+                # a c/(x - oo)^j term in ybar for some c and j
+                if poles[i] == oo and c[i][choice[i + 1]][j] != 0:
+                    return m, ybar, False
+                ybar += c[i][choice[i + 1]][j]/(x - poles[i])**(j+1)
+            m -= c[i][choice[i + 1]][0]
+
+        # Calculate the second summation for ybar
+        for i in range(N+1):
+            ybar += d[choice[0]][i]*x**i
+        return m, ybar, True
+
+    def construct_d(self, a, x, val_inf):
+        ser = a.series(x, oo)
+        N = -val_inf//2
+
+        # Case 4
+        if val_inf < 0:
+            temp = [0 for i in range(N+2)]
+            temp[N] = sqrt(ser.coeff(x, 2*N))
+            for s in range(N-1, -2, -1):
+                sm = 0
+                for j in range(s+1, N):
+                    sm += temp[j]*temp[N+s-j]
+                if s != -1:
+                    temp[s] = (ser.coeff(x, N+s) - sm)/(2*temp[N])
+            temp1 = list(map(lambda x: -x, temp))
+            temp[-1] = (ser.coeff(x, N+s) - temp[N] - sm)/(2*temp[N])
+            temp1[-1] = (ser.coeff(x, N+s) - temp1[N] - sm)/(2*temp1[N])
+            d = [temp, temp1]
+
+        # Case 5
+        elif val_inf == 0:
+            # List to store coefficients for plus case
+            temp = [0, 0]
+
+            # d_0  = sqrt(a_0)
+            temp[0] = sqrt(ser.coeff(x, 0))
+
+            # d_(-1) = a_(-1)/(2*d_0)
+            temp[-1] = ser.coeff(x, -1)/(2*temp[0])
+
+            # Coefficients for the minus case are just the negative
+            # of the coefficients for the positive case. Append both.
+            d = [temp, list(map(lambda x: -x, temp))]
+
+        # Case 6
+        else:
+            # s_oo = lim x->0 1/x**2 * a(1/x) which is equivalent to
+            # s_oo = lim x->oo x**2 * a(x)
+            s_inf = limit(x**2*a, x, oo)
+
+            # d_(-1) = (1 +- sqrt(1 + 4*s_oo))/2
+            d = [[(1 + sqrt(1 + 4*s_inf))/2], [(1 - sqrt(1 + 4*s_inf))/2]]
+        return d
+
+    def construct_c(self, a, x, poles, muls):
+        c = []
+        for pole, mul in zip(poles, muls):
+            c.append([])
+
+            # Case 3
+            if mul == 1:
+                # If multiplicity is 1, the coefficient to be added
+                # in the c-vector is 1 (no choice)
+                c[-1].extend([[1], [1]])
+
+            # Case 1
+            elif mul == 2:
+                # Find the coefficient of 1/(x - pole)**2 in the
+                # Laurent series expansion of a(x) about pole.
+                r = a.series(x, pole).coeff(x - pole, -2)
+
+                # If multiplicity is 2, the coefficient to be added
+                # in the c-vector is c = (1 +- sqrt(1 + 4*r))/2
+                c[-1].extend([[(1 + sqrt(1 + 4*r))/2], [(1 - sqrt(1 + 4*r))/2]])
+
+            # Case 2
+            else:
+                # Generate the coefficients using the recurrence
+                # relation mentioned in (5.14) in the thesis (Pg 80)
+
+                # r_i = mul/2
+                ri = mul//2
+
+                # Find the Laurent series about the pole
+                ser = a.series(x, pole)
+
+                # Start with an empty memo to store the coefficients
+                # This is for the plus case
+                temp = [0 for i in range(ri)]
+
+                # Base Case
+                temp[ri-1] = sqrt(ser.coeff(x - pole, -mul))
+
+                # Iterate backwards to find all coefficients
+                for s in range(ri-1, 0, -1):
+                    sm = 0
+                    for j in range(s+1, ri):
+                        sm += temp[j-1]*temp[ri+s-j-1]
+                    if s!= 1:
+                        temp[s-1] = (ser.coeff(x - pole, -ri-s) - sm)/(2*temp[ri-1])
+
+                # Memo for the minus case
+                temp1 = list(map(lambda x: -x, temp))
+
+                # Find the 0th coefficient in the recurrence
+                temp[0] = (ser.coeff(x - pole, -ri-s) - sm + ri*temp[ri-1])/(2*temp[ri-1])
+                temp1[0] = (ser.coeff(x - pole, -ri-s) - sm + ri*temp1[ri-1])/(2*temp1[ri-1])
+
+                # Add both the plus and minus cases' coefficients
+                c[-1].extend([temp, temp1])
+        return c
+
+    def find_poles(self, a, x):
+        # All roots of denominator are poles of a(x)
+        p = roots(denom(a), x)
+
+        # Substitute 1/x for x and check if oo is a pole
+        a = cancel(a.subs(x, 1/x).together())
+
+        # Find the poles of a(1/x)
+        p_inv = roots(denom(a))
+
+        # If 0 is a pole of a(1/x), oo is a pole of a(x)
+        if 0 in p_inv:
+            p.update({oo: p_inv[0]})
+        return p
+
+    def val_at_inf(self, a, x):
+        numer, denom = a.as_numer_denom()
+        # Valuation of a rational function at oo = deg(denom) - deg(numer)
+        return degree(denom, x) - degree(numer, x)
+
+    def _get_general_solution(self, *, simplify: bool = True):
+        # Step 0 :Match the equation
+        b0, b1, b2 = self.wilds_match()
+        fx = self.ode_problem.func
+        x = self.ode_problem.sym
+
+        # Step 1 : Convert to Normal Form
+        a = -b0*b2 + b1**2/4 - b1.diff(x)/2 + 3*b2.diff(x)**2/(4*b2**2) + b1*b2.diff(x)/(2*b2) - b2.diff(x, 2)/(2*b2)
+        a_t = cancel(a.together())
+
+        # Step 2
+        presol = []
+
+        # Step 3 : a(x) is 0
+        if a_t == 0:
+            presol.append(1/(x + self.ode_problem.get_numbered_constants(num=1)[0]))
+
+        # Step 4 : a(x) is a non-zero constant
+        elif x not in a_t.free_symbols:
+            presol.extend([sqrt(a), -sqrt(a)])
+
+        # Step 5 : Find poles and valuation at infinity
+        poles = self.find_poles(a_t, x)
+        poles, muls = list(poles.keys()), list(poles.values())
+        val_inf = self.val_at_inf(a_t, x)
+
+        if len(poles):
+            # Check necessary conditions (outlined in the module docstring)
+            if val_inf%2 != 0 or not all(map(lambda mul: (mul == 1 or (mul%2 == 0 and mul >= 2)), muls)):
+                raise ValueError("Rational Solution doesn't exist")
+
+            # Step 6
+            # Construct c-vectors for each singular point
+            c = self.construct_c(a, x, poles, muls)
+
+            # Construct d vectors for each singular point
+            d = self.construct_d(a, x, val_inf)
+
+            # Step 7 : Iterate over all possible combinations and return solutions
+            for it in range(2**(len(poles) + 1)):
+                # For each possible combination, generate an array of 0's and 1's
+                # where 0 means pick 1st choice and 1 means pick the second choice.
+                choice = list(map(lambda x: int(x), bin(it)[2:].zfill(len(poles) + 1)))
+
+                # Step 8 and 9 : Compute m and ybar
+                m, ybar, exists = self.compute_degree(x, poles, choice, c, d, -val_inf//2)
+
+                # Step 10 : Check if a valid solution exists. If yes, also check
+                # if m is a non-negative integer
+                if exists and S(m).is_nonnegative == True and S(m).is_integer == True:
+
+                    # Step 11 : Find polynomial solutions of degree m for the auxiliary equation
+                    psol, coeffs, exists = self.solve_aux_eq(a, x, m, ybar)
+
+                    # Step 12 : If valid polynomial solution exists, append solution.
+                    if exists:
+                        # m == 0 case
+                        if psol == 1 and coeffs == 0:
+                            # p(x) = 1, so p'(x)/p(x) term need not be added
+                            presol.append(ybar)
+                        # m is a positive integer and there are valid coefficients
+                        elif len(coeffs):
+                            # Substitute the valid coefficients to get p(x)
+                            psol = psol.subs(coeffs[0])
+                            # y(x) = ybar(x) + p'(x)/p(x)
+                            presol.append(ybar + psol.diff(x)/psol)
+
+                # If m is an expression with symbols, we generate a Piecewise solution where
+                # the solution will exist only if the expression is a nonnegative integer.
+                # NOTE: In this case, it seems like x**m + C1 is always a solution
+                # for the auxiliary equation. It is however NOT mentioned in the thesis.
+                elif len(m.free_symbols):
+                    C1 = Symbol('C1')
+                    psol = x**m + C1
+                    presol.append(ybar + psol.diff(x)/psol)
+        # Step 15 : Inverse transform the solutions of the equation in normal form
+        sol = list(map(lambda y: Eq(fx, -y/b2 - b2.diff(x)/(2*b2**2) - b1/(2*b2)), presol))
+        return sol
+
+# Avoid circular import:
+from .ode import get_numbered_constants
diff --git a/sympy/solvers/ode/tests/test_single.py b/sympy/solvers/ode/tests/test_single.py
index 8dd5f00c28d0..127c3a9e2dbb 100644
--- a/sympy/solvers/ode/tests/test_single.py
+++ b/sympy/solvers/ode/tests/test_single.py
@@ -528,6 +528,10 @@ def test_Riccati_special_minus2():
     _ode_solver_test(_get_examples_ode_sol_riccati)
 
 
+def test_1st_rational_riccati():
+    _ode_solver_test(_get_examples_ode_sol_1st_rational_riccati)
+
+
 def test_Bernoulli():
     _ode_solver_test(_get_examples_ode_sol_bernoulli)
 
@@ -780,6 +784,115 @@ def _get_examples_ode_sol_riccati():
     }
 
 
+@_add_example_keys
+def _get_examples_ode_sol_1st_rational_riccati():
+    a, b, c, k = symbols('a b c k')
+    # Type: 1st Order Rational Riccati, dy/dx = a + b*y + c*y**2,
+    # a, b, c are rational functions of x
+    return {
+            'hint': "1st_rational_riccati",
+            'func': f(x),
+            'examples':{
+    "rational_riccati_01": {
+        "eq": x**2 - (2*x + 1/x)*f(x) + f(x)**2 + Derivative(f(x), x),
+        "sol": [Eq(f(x), 2*x/(C1 + x**2) + (2*x**2 + 1)/(2*x) - 1/(2*x)), Eq(f(x), x)]
+    },
+    "rational_riccati_02": {
+        "eq": f(x)**2 + Derivative(f(x), x) - (4*x**6 - 8*x**5 + 12*x**4 + 4*x**3 + 7*x**2 - 20*x + 4)/(4*x**4),
+        "sol": [Eq(f(x), x + 2*x/(x**2 - 1) - 1 - 3/(2*x) + x**(-2))]
+    },
+    "rational_riccati_03": {
+        "eq": -x*f(x)**2 + Derivative(f(x), x) - 2*f(x)/x,
+        "sol": [Eq(f(x), -4*x**2/(C1 + x**4)), Eq(f(x), 0)]
+    },
+    "rational_riccati_04": {
+        "eq": -f(x)**2 + Derivative(f(x), x) + (15*x**2 - 20*x + 7)/((x - 1)**2*(2*x - 1)**2),
+        "sol": [Eq(f(x), Mul(-2, (2*x**3 - 5*x**2 + 4*x - 1), evaluate=False)/(C1 + x**4 - 10*x**3/3 + 4*x**2 - 2*x) +
+            (S(1)/2)*(1/(x - S(1)/2)) + 1/(x - 1)), Eq(f(x), (S(1)/2)*(1/(x - S(1)/2)) + 1/(x - 1))]
+    },
+    "rational_riccati_05": {
+        "eq": 9*x**2/4 - f(x)**2 + Derivative(f(x), x) - S(21)/2,
+        "sol": [Eq(f(x), 3*x/2 - (3*x**2 - 1)/(x*(x**2 - 1)))]
+    },
+    "rational_riccati_06": {
+        "eq": Derivative(f(x), x) + (3*x**2 + 1)*f(x)**2/x + (6*x**2 - x + 3)*f(x)/(x*(x - 1)) + (3*x**2 - 2*x + 2)/(x*(x - 1)**2),
+        "sol": [Eq(f(x), (-x**2*(-6 + (3*x**2 + 1)/x**2)/Mul(2, (3*x**2 + 1), evaluate=False) - x*(-3/Mul(2, (x + sqrt(3)*I/3), \
+            evaluate=False) + 1/Mul(2, (x - sqrt(3)*I/3), evaluate=False) - 1/(x - 2*sqrt(3)*I/3)) - (6*x**2 - x + 3)/Mul(2, (x - 1), \
+            evaluate=False))/(3*x**2 + 1)), Eq(f(x), (-x**2*(-6 + (3*x**2 + 1)/x**2)/Mul(2, (3*x**2 + 1), evaluate=False) - \
+            x*(-1/(x + 2*sqrt(3)*I/3) + 1/Mul(2, (x + sqrt(3)*I/3), evaluate=False) - 3/Mul(2, (x - sqrt(3)*I/3), evaluate=False)) \
+            - (6*x**2 - x + 3)/Mul(2, (x - 1), evaluate=False))/(3*x**2 + 1)), Eq(f(x), (-x**2*(-6 + (3*x**2 + 1)/x**2)/Mul(2, \
+            (3*x**2 + 1), evaluate=False) - x*(Mul(-1, (3*x**2 + 1), evaluate=False)/(C1 + x**3 + x) + 1/Mul(2, (x + sqrt(3)*I/3), \
+            evaluate=False) + 1/Mul(2, (x - sqrt(3)*I/3), evaluate=False)) - (6*x**2 - x + 3)/Mul(2, (x - 1), evaluate=False)) \
+            /(3*x**2 + 1)), Eq(f(x), (-x**2*(-6 + (3*x**2 + 1)/x**2)/(3*x**2 + 1) - x*(1/(x + sqrt(3)*I/3) + 1/(x - sqrt(3)*I/3)) \
+            - (6*x**2 - x + 3)/(x - 1))/Mul(2, (3*x**2 + 1), evaluate=False))]
+    },
+    "rational_riccati_07": {
+        "eq": f(x)**2 + Derivative(f(x), x) - 15/(4*x**2),
+        "sol": [Eq(f(x), 4*x**3/(C1 + x**4) - 3/(2*x)), Eq(f(x), -3/(2*x))]
+    },
+    "rational_riccati_08": {
+        "eq": 3*f(x)**2 + Derivative(f(x), x) - 2/x**2,
+        "sol": [Eq(f(x), 5*x**4/Mul(3, (C1 + x**5), evaluate=False) - 2/(3*x)), Eq(f(x), -2/(3*x))]
+    },
+    "rational_riccati_09": {
+        "eq": -a*b*k**2*f(x)**2 + 2*a*b*k*f(x) - a*b + Derivative(f(x), x),
+        "sol": [Eq(f(x), (1 - 1/(a*b*k*(C1 + x)))/k)]
+    },
+    "rational_riccati_10": {
+        "eq": Derivative(f(x), x) - 2*I*(f(x)**2 + 1)/x,
+        "sol": [Eq(f(x), -I*x*(-4*x**3/(C1 + x**4) + 3/(2*x))/2 - I/4), Eq(f(x), -I)]
+    },
+    "rational_riccati_11": {
+        "eq": Derivative(f(x), x) - f(x)**2/x + 1/x,
+        "sol": [Eq(f(x), x*(-2*x/(C1 + x**2) + 1/(2*x)) + S(1)/2), Eq(f(x), 1)]
+    },
+    "rational_riccati_12": {
+        "eq": a*f(x)**2 + c + x*Derivative(f(x), x),
+        "sol": [Eq(f(x), sqrt(-a*c)*(2*x**(2*sqrt(-a*c))/(C1 + x**(2*sqrt(-a*c))) - 1)/a), Eq(f(x), -sqrt(-a*c)/a)]
+    },
+    "rational_riccati_13": {
+        "eq": a*f(x)**2/x - b*f(x)/x - c/x + Derivative(f(x), x),
+        "sol": [Eq(f(x), (b/2 + x**(sqrt(4*a*c + b**2))*sqrt(4*a*c + b**2)/(C1 + x**(sqrt(4*a*c + b**2))) - \
+            sqrt(4*a*c + b**2)/2)/a), Eq(f(x), (b - sqrt(4*a*c + b**2))/(2*a))]
+    },
+    "rational_riccati_14": {
+        "eq": f(x)**2 + Derivative(f(x), x) - f(x)/x,
+        "sol": [Eq(f(x), 2*x/(C1 + x**2)), Eq(f(x), 0)]
+    },
+    "rational_riccati_15": {
+        "eq": -x**2 - (2*x + 1/x)*f(x) - f(x)**2 + Derivative(f(x), x),
+        "sol": [Eq(f(x), -2*x/(C1 + x**2) - (2*x**2 + 1)/(2*x) + 1/(2*x)), Eq(f(x), -x)]
+    },
+    "rational_riccati_16": {
+        "eq": f(x)**2 + Derivative(f(x), x) + 4*f(x)/x + 2/x**2,
+        "sol": [Eq(f(x), 1/(C1 + x) - 2/x)]
+    },
+    "rational_riccati_17": {
+        "eq": a*f(x)**2 - b/x**2 + Derivative(f(x), x),
+        "sol": [Eq(f(x), Mul(-1, (-x**(sqrt(4*a*b + 1))*sqrt(4*a*b + 1)/(C1 + x**(sqrt(4*a*b + 1))) + \
+            sqrt(4*a*b + 1)/2 - S(1)/2), evaluate=False)/(a*x)), Eq(f(x), (1 - sqrt(4*a*b + 1))/(2*a*x))]
+    },
+    "rational_riccati_18": {
+        "eq": a*f(x)**2 + Derivative(f(x), x) + (b + c)/x**2,
+        "sol": [Eq(f(x), Mul(-1, (-x**(sqrt(-4*a*(b + c) + 1)/2 + sqrt(-4*a*b - 4*a*c + 1)/2)*(sqrt(-4*a*(b + c) + 1) \
+            + sqrt(-4*a*b - 4*a*c + 1))/(C1 + x**(sqrt(-4*a*(b + c) + 1)/2 + sqrt(-4*a*b - 4*a*c + 1)/2)) + sqrt(-4*a*b \
+            - 4*a*c + 1) - 1), evaluate=False)/(2*a*x)), Eq(f(x), Mul(-1, (-x**(-sqrt(-4*a*(b + c) + 1)/2 + sqrt(-4*a*b \
+            - 4*a*c + 1)/2)*(-sqrt(-4*a*(b + c) + 1) + sqrt(-4*a*b - 4*a*c + 1))/(C1 + x**(-sqrt(-4*a*(b + c) + 1)/2 + \
+            sqrt(-4*a*b - 4*a*c + 1)/2)) + sqrt(-4*a*b - 4*a*c + 1) - 1), evaluate=False)/(2*a*x))]
+    },
+    "rational_riccati_19": {
+        "eq": 2*x**2*Derivative(f(x), x) - x*(4*f(x) + Derivative(f(x), x) - 4) + (f(x) - 1)*f(x),
+        "sol": [Eq(f(x), -x*(1 - 2*x)/(C1 + x) - 2*x - (S(1)/2 - x)*(4*x + 1)/(2*x - 1) + S(1)/2)]
+    },
+    "rational_riccati_20": {
+        "eq": x**4*Derivative(f(x), x) + x**2 - x*(2*f(x)**2 + Derivative(f(x), x)) + f(x),
+        "sol": [Eq(f(x), 3*x**2/4 + (x**3 - 1)*(-2*x/(C1 + x**2) + 1/(2*x))/2 + 1/(4*x)), \
+            Eq(f(x), 3*x**2/4 + (x**3 - 1)/(4*x) + 1/(4*x))]
+    }
+    }
+    }
+
+
 @_add_example_keys
 def _get_examples_ode_sol_1st_linear():
     # Type: first order linear form f'(x)+p(x)f(x)=q(x)


Code review comments:
diff_hunk: @@ -0,0 +1,500 @@
+r"""
+This module contains :py:meth:`~sympy.solvers.ode.riccati.RationalRiccati`,
+a solver class which gives rational solutions to first order Riccati ODEs.
+A general first order Riccati ODE is given by
+
+.. math:: y' = b_0(x) + b_1(x)w + b_2(x)w^2
+
+where `b_0, b_1` and `b_2` can be arbitrary functions of `x` with `b_2 \ne 0`.
+
+Background
+==========
+
+A Riccati equation can be transformed to its normal form
+
+.. math:: y' + y^2 = a(x)
+
+using the transformation
+
+.. math:: y = -b_2(x) - \frac{b'_2(x)}{2 b_2(x)} - \frac{b_1(x)}{2}
+
+where `a(x)` is given by
+
+.. math:: a(x) = \frac{1}{4}\left(\frac{b_2'}{b_2} + b_1\right)^2 - \frac{1}{2}\left(\frac{b_2'}{b_2} + b_1\right)' - b_0 b_2
+
+Thus, we can develop an algorithm to solve for the Riccati equation
+in its normal form, which would in turn give us the solution for
+the original Riccati equation.
+
+Algorithm
+=========
+
+The algorithm implemented here was presented in a Ph.D thesis by
+N. Thieu Vo. The entire thesis can be found here -
+https://www3.risc.jku.at/publications/download/risc_5387/PhDThesisThieu.pdf
+
+We have only implemented the Rational Riccati solver (Pg 78 in Thesis).
+Before we proceed towards the implementation of the algorithm, a
+few definitions to understand -
+
+1. Valuation of a Rational Function at `\infty`:
+    The valuation of a rational function `p(x)` at `\infty` is equal
+    to the difference between the degree of the denominator and the
+    numerator of `p(x)`.
+
+    Note: A general definition of valuation of a rational function
+    at any value of `x` can be found in Pg 63 of the thesis, but
+    is not of any interest for this algorithm.
+
+2. Zeros and Poles of a Rational Function:
+    Let `a(x) = \frac{S(x)}{T(x)}, T \ne 0` be a rational function
+    of `x`. Then -
+
+    a. The Zeros of `a(x)` are the roots of `S(x)`.
+    b. The Poles of `a(x)` are the roots of `T(x)`. However, `\infty`
+    can also be a pole of a(x). We say that `a(x)` has a pole at
+    `\infty` if `a(\frac{1}{x})` has a pole at 0.
+
+Every pole is associated with an order that is equal to the multiplicity
+of its appearence as a root of `T(x)`. A pole is called a simple pole if
+it has an order 1. Similarly, a pole is called a multiple pole if it has
+an order `\ge` 2.
+
+Necessary Conditions
+====================
+
+For a Riccati equation in its normal form,
+
+.. math:: y' + y^2 = a(x)
+
+we can define
+
+a. A pole is called a movable pole if it is a pole of `y(x)` and is not
+a pole of `a(x)`.
+b. Similarly, a pole is called a non-movable pole if it is a pole of both
+`y(x)` and `a(x)`.
+
+Then, the algorithm states that a rational solution exists only if -
+
+a. Every pole of `a(x)` must be either a simple pole or a multiple pole
+of even order.
+b. The valuation of `a(x)` at `\infty` must be even or be `\ge` 2.
+
+Implementation
+==============
+
+The code is written to match the steps given in the thesis (Pg 82)
+
+Step 0 : Match the equation -
+Find `b_0, b_1` and `b_2`. If `b_2 = 0` or no such functions exist, raise
+an error
+
+Step 1 : Transform the equation to its normal form as explained in the
+theory section.
+
+Step 2 : Initialize 2 empty set of solutions, ``sol`` and ``presol``.
+
+Step 3 : If `a(x) = 0`, append `\frac{1}/{(x - C1)}` to ``presol``.
+
+Step 4 : If `a(x)` is a rational non-zero number, append `\pm \sqrt{a}`
+to ``presol``.
+
+Step 5 : Find the poles and their multiplicities of `a(x)` using
+``find_poles``. Let the number of poles be `n`. Also find the
+valuation of `a(x)` at `\infty` using ``val_at_inf``.
+
+Step 6 : Find `n` c-vectors (one for each pole) and 1 d-vector using
+``construct_c`` and ``construct_d``. Now, determine all the ``2**(n + 1)``
+combinations of choosing between 2 choices for each of the `n` c-vectors
+and 1 d-vector.
+
+For each of these above combinations, do
+
+Step 8 : Compute `m` in ``compute_degree``
+
+Step 9 : In ``compute_degree``, compute ybar as well. If ybar contains any
+term like `\frac{c_{ij}}{(x - oo)^j}` where `c_{ij}` is non-zero, discard
+this solution.
+
+Step 10 : If `m` is a non-negative integer -
+
+Step 11: Find a polynomial solution of degree `m` for the auxiliary equation.
+
+There are 2 cases possible -
+
+    a. `m` is a non-negative integer: We can solve for the coefficients
+    in `p(x)` using Undetermined Coefficients.
+
+    b. `m` is an expression: In this case, we may not be able to determine
+    if it is possible for `m` to be a non-negative integer. However, the
+    solution to the auxiliary equation always seems to be x**m + C1 (where
+    C1 is an arbitrary constant). This has been confirmed by testing many
+    examples. NOTE that this is NOT mentioned in the thesis and is an
+    assumption made by me. If any failing test cases are found, this case
+    must be removed.
+
+Step 12 : For each `p(x)` that exists, append `ybar + \frac{p'(x)}{p(x)}`
+to ``presol``.
+
+Step 13 : For each solution in ``presol``, apply an inverse transformation
+and append them to ``sol``, so that the solutions of the original equation
+are found using the solutions of the equation in its normal form.
+"""
+
+from .single import SinglePatternODESolver
+
+from sympy.core import S
+from sympy.core.add import Add
+from sympy.core.function import Function
+from sympy.core.numbers import oo
+from sympy.core.relational import Eq
+from sympy.core.symbol import Symbol, Wild
+from sympy.functions import sqrt
+from sympy.polys.polytools import cancel, degree
+from sympy.polys.polyroots import roots
+from sympy.simplify.radsimp import numer, denom
+from sympy.series.limits import limit
+from sympy.solvers.solvers import solve
+
+class RationalRiccati(SinglePatternODESolver):
path: sympy/solvers/ode/riccati.py
text: @user1:
I don't think this should use SinglePatternODESolver.

Rather the code here should just be a bunch of functions for doing the different parts of the algorithm. Then we can make a Solve class in single.py that wraps this for the current dsolve framework.

@author:
So we don't need the `RationalRiccati` class? What would the `Solve` class constitute?

@user1:
The solver class would basically be this one except that it wouldn't have all of these methods. Instead it would just call the functions from the ricatti module.

These classes exist just for the current pattern-matching structure of dsolve but we might want to move beyond that to a different structure later so an algorithm like this should not be tied to that structure.