Calculating the transfer function from electronic and mechanical equations in sympy

I am trying to calculate the transfer function of an electric motor using sympy in an effort to save me time doing hand calculations. Here are the mechanical and electrical equations for the motor:

Motor equations

Motor Electrical Equation: V(t) - L di/dt - i(t)R - e(t) = 0

Motor Mechanical Equation: J * d^2theta/dt^2 = T_m - b * omega

Motor Torque T_m = k * i

I have then turned these into Sympy equations

import sympy as sym

v_in, L, R, e_back, i, t, j, theta, damp, T_s , omega, T_m, k = 
sym.symbols(r"V_in, L, R, e_back, i, t, J, theta, b, T_shaft, \dot{\theta}, T_motor, k", real=True)
v_out = L * Derivative(i, t) - i * R - e_back
i_t = Derivative(i, t)

motor_elec_eqn = sym.Eq(v_in, v_out)
display(motor_elec_eqn)

theta_dot = Derivative(j, t)
theta_dot_dot = Derivative(theta_dot, t)

torque_shaft = sym.Eq(j * theta_dot_dot, T_s - damp * theta_dot)
torque_motor = sym.Eq(T_s,k * i)

tf = sym.solve(theta_dot,theta_dot_dot,torque_motor,torque_shaft,motor_elec_eqn)
sym.pretty(tf)

This method gives me an empty tf, i.e. tf = []

I then tried to convert the original equations into the Laplace domain in the hope that might produce something that could be solved:

import sympy as sym

# Define symbolic variables
s, V_in, L, R, e_back, i, J, theta, b, T_s, omega, T_m, k = sym.symbols(
    r"s, V_in, L, R, e_back, i, J, theta, b, T_s, omega, T_m, k", real=True
)

# Define Laplace transforms of variables
I_s = sym.LaplaceTransform(i, t, s)
Theta_s = sym.LaplaceTransform(theta, t, s)
E_back_s = sym.LaplaceTransform(e_back, t, s)
T_s_s = sym.LaplaceTransform(T_s, t, s)
V_out_s = L * (s * I_s) - R * I_s - E_back_s
display(V_out_s)
Motor_elec_eqn_s = sym.Eq(V_in, V_out_s)

# Define Laplace domain equations for the mechanical system
Theta_dot_s = s * Theta_s
Theta_dot_dot_s = s * Theta_dot_s
Torque_shaft_s = sym.Eq(J * Theta_dot_dot_s,  -b * Theta_dot_s + T_s_s)
Torque_motor_s = sym.Eq(T_m, k * I_s)

display(Motor_elec_eqn_s)
display(Torque_motor_s)
display(Torque_shaft_s)

# Solve for the transfer function
transfer_function = sym.solve(
    [Motor_elec_eqn_s, Torque_shaft_s, Torque_motor_s], Theta_s
)[0]

# Simplify the transfer function if needed
transfer_function = sym.simplify(transfer_function)

# Display the transfer function
sym.pretty(transfer_function)

That gave a different error that I don't understand and don't know how to fix (See image). Is it possible to get the transfer function using sympy, any other method or will I need to hand calculate the results and input those into sympy?

I've noticed a few errors in your equations, so I'm going to use a slightly different notation because I'm following step by step the procedure shown on the book in my possession.

I'm also going to use Algebra with SymPy, a package that exposes a very handy Equation class which can be used with math operators (+-*/).

from sympy import *
from algebra_with_sympy import Equation, solve, algwsym_config
algwsym_config.output.label = False
algwsym_config.output.solve_to_list = True

Ra, La, Jm, Dm, Kb, Kt = symbols("R_a, L_a, J_m, D_m, K_b, K_t")
s, t = symbols("s, t")
ia, Ea, Vb, Tm, theta_m = [symbols(n, cls=Function)(t) for n in ["i_a", "E_a", "V_b", "T_m", "theta_m"]]

# back-electromotive force
vb_eq = Equation(Vb, Kb * theta_m.diff(t))
# circuit equation
ci_eq = Equation(Ra * ia + La * ia.diff(t) + Vb, Ea)
# torque equation
torque_eq = Equation(Tm, Kt * ia)
# mechanical equation
me_eq = Equation(Tm, Jm * theta_m.diff(t, 2) + Dm * theta_m.diff(t))

display(vb_eq, ci_eq, torque_eq, me_eq)

Let's convert these time-domain equation to the Laplace domain.

def laplace_transform_eq(eq):
    # laplace_transform is not able to deal with objects of type `Equation`
    # or `Eq`, so we have to manually apply it to both sides of the equation.
    new_eq = Equation(*[laplace_transform(term.doit(), t, s, noconds=True) for term in eq.args])
    # apply zero initial conditions
    functions_of_t = eq.find(Function)
    d = {}
    for f in functions_of_t:
        d.update({f.subs(t, 0): 0, f.diff(t).subs(t, 0): 0})
    return new_eq.subs(d)

vb_eq = laplace_transform_eq(vb_eq)
ci_eq = laplace_transform_eq(ci_eq)
torque_eq = laplace_transform_eq(torque_eq)
me_eq = laplace_transform_eq(me_eq)

display(vb_eq, ci_eq, torque_eq, me_eq)

The above results really hurts my eyes, so I'm going to replace the Laplace transforms with ordinary symbols. Which leads me to this recommendation: if possible, write your equations in the Laplace domain, so you save yourself some typing:

laplace_terms = [laplace_transform(f, t, s, noconds=True) for f in [Vb, Tm, ia, Ea, theta_m]]
Vb, Tm, ia, Ea, theta_m = symbols("V_b, T_m, i_a, E_a, theta_m")
d = {k: v for k, v in zip(laplace_terms, [Vb, Tm, ia, Ea, theta_m])}
vb_eq, ci_eq, torque_eq, me_eq = [eq.subs(d) for eq in [vb_eq, ci_eq, torque_eq, me_eq]]

display(vb_eq, ci_eq, torque_eq, me_eq)

Now, we substitutes some equation into some other equation. I'm sure you have the manual procedure listed in your book:

ia_eq = solve(torque_eq, ia)[0]
ia_eq

tf = (ci_eq.subs(ia_eq, vb_eq).collect(Tm / Kt)).subs(me_eq).collect(theta_m)
tf

tf = (1 / (tf / theta_m)).swap
tf

Finally, I extract the right-hand-side of the transfer function and apply some other manipulations:

tf.rhs.together().expand()