3.41 function to solve reduced Riccati ode in Maple?

This is a standalone function that solves \(y'=a x^n + b y^2\). This is called the specific or reduced Riccati ode. The formulas used are from Eqworld ode0106 and Dr Dobrushkin web page

#solves y' = a*x^n + b*y^2, using standard formulas for n=-2 and n<>-2 
 
special_riccati_dsolve := proc(n,a,b,func::function(name),$) 
 
    local y:=op(0,func); 
    local x:=op(1,func); 
    local sol; 
    local C1,C2; 
 
    if not (n::'integer' or n::'fraction') then 
       print("argument n  must be integer or fraction"); 
       RETURN(NULL); 
    fi; 
 
    if n=-2 then 
        proc() 
            local z,sol_z,lambda; 
            try 
                sol_z:=timelimit(30,[solve(b*z^2+z+a=0,z)]); 
                if nops(sol_z)=0 then 
                   error "Unable to solve for root"; 
                fi; 
            catch: 
                error "Unable to solve for root"; 
            end try; 
 
            z:=sol_z[1]; #pick one, any will work 
            sol:=y(x)= z/x - x^(2*b*z)/( b*x/(2*b*z+1)*x^(2*b*z) + c__1); 
        end proc(); 
    else 
        proc() 
            local k,w; 
            k:=1+n/2; 
            C1:= c__1; 
            C2:= c__2; 
            try 
                if evalb(a*b >0) then 
                    w:=  sqrt(x) * ( C1*BesselJ(1/(2*k),1/k*sqrt(a*b)*x^k) + C2* BesselY(1/(2*k),1/k*sqrt(a*b)*x^k) ); 
                else 
                    w:=  sqrt(x) * ( C1*BesselI(1/(2*k),1/k*sqrt(-a*b)*x^k) + C2* BesselK(1/(2*k),1/k*sqrt(-a*b)*x^k) ); 
                fi; 
            catch: 
                w:=  sqrt(x) * ( C1*BesselJ(1/(2*k),1/k*sqrt(a*b)*x^k) + C2* BesselY(1/(2*k),1/k*sqrt(a*b)*x^k) ); 
            end try; 
            sol:= simplify(-1/b*diff(w,x)/w); 
            sol:= y(x)=eval(sol,C2=1); 
        end proc(); 
    fi; 
 
    RETURN(simplify(sol)); 
 
end proc:
 

Example usages

a:=1:n:=-2:b:=1: 
ode:=diff(y(x),x)=a*x^n+b*y(x)^2; 
sol:=special_riccati_dsolve(n,a,b,y(x)); 
 
    sol := y(x) = ((-x^(sqrt(3)*I) - 3*c__1)*sqrt(3)*I + 3*x^(sqrt(3)*I) + 3*c__1)/(2*(x^(sqrt(3)*I)*sqrt(3)*I - 3*c__1)*x) 
 
a:=5:n:=-8/3:b:=3: 
ode:=diff(y(x),x)=a*x^n+b*y(x)^2; 
sol:=special_riccati_dsolve(n,a,b,y(x)); 
 
   sol := y(x) = -sqrt(15)*((sqrt(15)*x - 15*c__1*x^(2/3) + c__1*x^(4/3)/3)*cos(3*sqrt(15)/x^(1/3)) + sin(3*sqrt(15)/x^(1/3))*(sqrt(15)*c__1*x + 15*x^(2/3) - x^(4/3)/3))/(((c__1*x^(1/3)*sqrt(15) + 45)*cos(3*sqrt(15)/x^(1/3)) + 45*(-sqrt(15)*x^(1/3)/45 + c__1)*sin(3*sqrt(15)/x^(1/3)))*x^2) 
 
 
 
a:=5:n:=5:b:=2: 
ode:=diff(y(x),x)=a*x^n+b*y(x)^2; 
sol:=special_riccati_dsolve(n,a,b,y(x)); 
 
   sol := y(x) = -sqrt(10)*x^(5/2)*(BesselY(-6/7, (2*sqrt(10)*x^(7/2))/7) + BesselJ(-6/7, (2*sqrt(10)*x^(7/2))/7)*c__1)/(2*c__1*BesselJ(1/7, (2*sqrt(10)*x^(7/2))/7) + 2*BesselY(1/7, (2*sqrt(10)*x^(7/2))/7))
 

Ofcourse, dsolve can be used to solve these also.