Generate fractions in Taylor series coefficients a bit more robustly

Author: hersle

docs/src/tutorials/perturbation.md

@@ -79,8 +79,8 @@ function solve_perturbed(eq, x, x₀, ϵ, order)
 
     # solve higher-order equations order-by-order
     for i in 2:length(eqs)
+        eqs[i] = substitute(eqs[i], sol) # substitute lower-order solutions
         x_coeff = Symbolics.symbolic_linear_solve(eqs[i], x_coeffs[i]) # solve linear n-th order equation for x_n
-        x_coeff = substitute(x_coeff, sol) # substitute lower-order solutions to get numerical value
         sol = merge(sol, Dict(x_coeffs[i] => x_coeff)) # store solution
     end
 

src/taylor.jl

@@ -72,7 +72,7 @@ function taylor_coeff(f, x, n = missing; rationalize=true)
     if rationalize && unwrap(c) isa Number
         # TODO: make rational coefficients "organically" and not using rationalize (see https://github.com/JuliaSymbolics/Symbolics.jl/issues/1299)
         c = unwrap(c)
-        c = isinteger(c) ? Int(c) : Base.rationalize(c) # convert non-integer floats to rational numbers; avoid name clash between rationalize and Base.rationalize()
+        c = Base.rationalize(c) # convert integers/floats to rational numbers; avoid name clash between rationalize and Base.rationalize()
     end
     return c
 end