Analytical expressions for ODE integration methods using Sympy

I am trying to use Sympy to calculate analytical expressions for the error of numerical solutions of ordinary differential equations (ODEs), such as the Euler and Runge-Kutta methods.

We have the ODE: $\dot x = f(x)$ with initial value $x(t_0) = x_0$, and we want to find the error of the approximate solution $x(t+\Delta t) \approx x_0 + \Delta t f(x_0) + \Delta t^2 f'(x_0)f(x_0)/2 +...$.

The procedure is to expand the solution $x(t)$ in a Taylor series up to two orders higher than the numerical method under analysis and then compare this expansion to the expanded evaluation of the numerical solution. For the Euler method this is very straightforward (easy), because both expansions don't require much algebra:

Euler method: x1 = x0 + \Delta t f(x_0)

Taylor series, up to $O(\Delta t^2)$: $x(t+\Delta t) = x_0 + \Delta t f(x_0) + \Delta t^2 f'(x_0)f(x_0)/2 + O(\Delta t^3) $

Error: $|E| = | \Delta t^2 f'(x_0)f(x_0)/2 | $.

The trick above is to use $dx/dt = f(x)$ and $d^2x/dt^2 = (df/dx)(dx/dt)$ (chain rule) in the Taylor expansion of $x(t)$. However, in order 2 or higher, things become more complicated, so I would like to use Sympy, to help in the development of the expressions.

Here is my attempt:

from sympy import *
init_printing()

t, t0, dt, x, x0, x1, x_dot = symbols("t t_0 {\\Delta}t x x_0 x_1 \\dot{x}") 
xt, f  = symbols("x f", cls=Function) #x(t) e f(x(t))
Erro =symbols("E")
fp = [f] # A Python list to contain the symbols of df/dx, up to order 9
string_fp = "f"
for i in range(1,10):
    string_fp +="'"
    fp.append(symbols(string_fp, cls=Function))

def my_taylor(f, x, x0, order):
    Taylor_expansion = f(x0)
    for i in range(1, order+1):
        Taylor_expansion += (1/factorial(i))*((x-x0)**i)* diff(f(x), x, i)
        #Taylor_expansion += (1/factorial(i))*((x-x0)**i)* diff(f(x), x, i).subs(x,x0)
        #Taylor_expansion += (1/factorial(i))*((x-x0)**i)*fp[i](x0)
    Taylor_expansion += Order((x-x0)**(order+1),(x,x0))
    return Taylor_expansion

#teste = my_taylor(f, x, x0, 2)

def my_other_taylor(xt, t, t0, f, order):
    Taylor_expansion = xt(t0)
    Taylor_expansion += (t-t0)*f(x0)
    for i in range(2, order+1):
        #Taylor_expansion += (1/factorial(i))*((t-t0)**i)* diff(f(x), x, i-1).subs(x,x0)
        #Taylor_expansion += (1/factorial(i))*((t-t0)**i)*fp[i](x0)
        Taylor_expansion += (1/factorial(i))*((t-t0)**i)* diff(f(xt(t)), t, i-1).subs(xt(t),x0)
    #Taylor_expansion += Order((factor(t-t0)**(order+1)))
    Taylor_expansion += Order((t-t0)**(order+1),(t,t0))
    return Taylor_expansion
#teste2 = my_other_taylor(xt, t, t0, f, 5)

#Taylor_f  = f(x).series(x, x0=x0, n=3)
#Taylor_x  = xt(t).series(t, x0=t0, n=3)
#Taylor_x_2 = Taylor_x.subs(diff(xt(t),t), f(x))  ## Esta substituição não está funcionando
Taylor_f = my_taylor(f, x, x0, 2)
Taylor_x = my_other_taylor(xt, t, t0, f, 4)
Taylor_x_dt = Taylor_x.subs(t, dt+t0)

#Taylor_x_2 = Taylor_x.subs(t-t0, dt)
pprint(Taylor_x_dt, use_unicode=True)

### Calcular x1 usando método de Euler
x1 = x0 + dt*f(x0)
Erro = abs(x1-Taylor_x_dt.removeO().subs(xt(t0),x0))

I am not getting the correct results, and I can see why, but I don't know how to fix the Taylor expansion. The derivatives of $f(x)$ are supposed to be "evaluated" at $x_0$, but I can't find how to make a symbolic evaluation of a derivative in sympy's documentation and tutorials. If I make .subs(x, x0), it substitutes the variable used in the differentiation , like $df(x)/dx$ becomes $df(x_0)/dx_0$, instead of $df(x_0)/dx$. Also, I can't find a way to apply the chain rule and substitute the time derivatives of $x(t)$ by the derivatives of $f$ w.r.t $x$.

source https://stackoverflow.com/questions/75881778/analytical-expressions-for-ode-integration-methods-using-sympy