How to use the sympy.simplify function in sympy

def add(p1, p2, s1=1, s2=1):
    x = p1[0] * s1 + p2[0] * s2
    y = p1[1] * s1 + p2[1] * s2
    return [simplify(x), simplify(y)]

print('%\\abs{B} =', Bmag)
    Wd_1 = Wd_1.subs(Bmag, 1/Binv)
    Wd_C = Wd_C.subs(Bmag, 1/Binv)
    Wd_S = Wd_S.subs(Bmag, 1/Binv)

    lhs = Wd_1 + Wd_C*BigC
    rhs = -Wd_S*BigS
    lhs = lhs**2
    rhs = rhs**2
    W = expand(lhs - rhs)
    W = expand(W.subs(1/Binv**2, Bmag**2))
    W = expand(W.subs(BigS**2, BigC**2 - 1))
    W = W.collect([BigC, BigC**2], evaluate=False)

    a = simplify(W[BigC**2])
    b = simplify(W[BigC])
    c = simplify(W[S.One])

    print('#%\\text{Require } aC^{2}+bC+c = 0')

    print('a =', a)
    print('b =', b)
    print('c =', c)

    x = Symbol('x')
    C = solve(a*x**2 + b*x + c, x)[0]
    print('%b^{2}-4ac =', simplify(b**2 - 4*a*c))
    print('%\\f{\\cosh}{\\alpha} = C = -b/(2a) =', expand(simplify(expand(C))))
    return

# create a small matrix M, and triangularize it
        M = create_ma(deg_f, deg_g, row1, row0, col_num)
        for i in range(deg_f - deg_g + 1):
            M1 = pivot(M, i, i)
            M = M1[:, :]

        # treat last row of M as poly; find its degree
        d = find_degree(M, deg_f)
        if d == None:
            return sr_list
        exp_deg = deg_g - 1

        # evaluate one determinant & make coefficients subresultants
        sign_value = correct_sign(n, m, s1, exp_deg, exp_deg - d)
        poly = row2poly(M[M.rows - 1, :], d, x)
        poly = simplify((poly / LC(poly, x)) * sign_value)

        # append poly with subresultant coeffs
        sr_list.append(poly)

        # update degrees and rows
        deg_f, deg_g = deg_g, d
        row0 = row1
        row1 = Poly(poly, x, domain = QQ).all_coeffs()
        leng1 = len(row1)
        for i in range(col_num - leng1):
            row1.append(0)
        row1 = Matrix([row1])

    return sr_list

alpha_star=sympy.Matrix(alpha[:-1,:])
        beta_star=sympy.Matrix(beta[:-1,:])
        I = sympy.eye(s)

        v_mp1 = v[-1]
        vstar = sympy.Matrix(v[:-1])
        alpha_mp1 = sympy.Matrix(alpha[-1,:]).T
        beta_mp1 = sympy.Matrix(beta[-1,:]).T

        if formula == 'det':
            xsym = I - alpha_star - z*beta_star + vstar/v_mp1 * (alpha_mp1+z*beta_mp1)
            p1 = sympy.simplify(xsym.det(method='berkowitz')*v_mp1)
            p1 = p1.as_poly(z).all_coeffs()

            denomsym = I - alpha_star - z*beta_star
            q1 = sympy.simplify(denomsym.det(method='berkowitz'))
            q1 = q1.as_poly(z).all_coeffs()

        elif formula == 'lts': # lower_triangular_solve
            p1 = (alpha_mp1 + z*beta_mp1)*((I-alpha_star-z*beta_star).lower_triangular_solve(vstar))
            p1 = sympy.poly(p1[0])+v_mp1
            p1 = p1.all_coeffs()

        elif formula == 'pow': # Power series
            apbz_star = alpha_star + beta_star*z
            apbz = sympy.Matrix(alpha_mp1+z*beta_mp1)

            # Compute (I-(alpha + z beta)^(-1) = I + (alpha + z beta) + (alpha + z beta)^2 + ... + (alpha + z beta)^(s-1)
            # This is coded for Shu-Osher coefficients
            # For them, we need to take m=s
            # For Butcher coefficients, perhaps we could get away with m=num_seq_dep_stages (???)
            apbz_power = I

### Compute normal to panel
# see surface.AeroGridSurface.get_panel_normal
R02=Zeta02-Zeta00
R13=Zeta03-Zeta01
Norm=linfunc.cross_product(R02,R13)
Norm=Norm/linfunc.norm2(Norm)
### check norm
assert linfunc.scalar_product(Norm,R02).simplify()==0, 'Norm is wrong'
assert linfunc.scalar_product(Norm,R13).simplify()==0, 'Norm is wrong'
assert linfunc.scalar_product(Norm,Norm).simplify()==1, 'Normal is not unit length'


### Compute normal velocity at panel
Unorm=linfunc.scalar_product(Norm,Uc)
Unorm=sm.simplify(Unorm)

### Compute derivative
dUnorm_dZeta=sm.derive_by_array(Unorm,[Zeta00,Zeta01,Zeta02,Zeta03])
#dUnorm_dZeta=linfunc.simplify(dUnorm_dZeta)



################################################################################
### exploit combined derivatives
################################################################################


dR_dZeta=sm.derive_by_array([R02,R13],[Zeta00,Zeta01,Zeta02,Zeta03])

### redefine R02,R13
r02_x,r02_y,r02_z=sm.symbols('r02_x r02_y r02_z', real=True)

# at point P so these internal forces will cancel. Note point P is fixed in
# both A and B.
forces = [(pO, K_EA), (pC_star, K_BC), (pB_hat, -K_BC),
          (pA_star, -mA*g*E.x), (pB_star, -mB*g*E.x),
          (pC_star, -mC*g*E.x), (pD_star, -mD*g*E.x)]
torques = [(A, T_EA - T_AB), (B, T_AB - T_BC), (C, T_BC)]

## --- define partial velocities ---
partials = partial_velocities([f[0] for f in forces + torques],
                              [u1, u2, u3], E, kde_map)

## -- calculate generalized active forces ---
Fr, _ = generalized_active_forces(partials, forces + torques)
print("Generalized active forces:")
for i, f in enumerate(Fr, 1):
    print("F{0} = {1}".format(i, msprint(simplify(f))))

# Don't pick a non-polynomial algebraic to be differentiated
            if rule == pull_out_algebraic and not u.is_polynomial(symbol):
                return
            # Don't trade one logarithm for another
            if isinstance(u, sympy.log):
                rec_dv = 1/dv
                if (rec_dv.is_polynomial(symbol) and
                    degree(rec_dv, symbol) == 1):
                        return

            for rule in liate_rules[index + 1:]:
                r = rule(integrand)
                # make sure dv is amenable to integration
                if r and r[0].subs(dummy, 1).equals(dv):
                    du = u.diff(symbol)
                    v_step = integral_steps(sympy.simplify(dv), symbol)
                    v = _manualintegrate(v_step)

                    return u, dv, v, du, v_step

def diff_eq(u, x):
        eqs =  {'diff': -sm.diff(u, x, x) - f,
                'BC1': sm.diff(u, x).subs(x, 0) - C,
                'BC2': u.subs(x, L) - D}
        for eq in eqs:
            eqs[eq] = sm.simplify(eqs[eq])

w = (E2|e3)
    w = w.expand()
    print('E2|e3 =',w)

    w = (E3|e1)
    w = w.expand()
    print('E3|e1 =',w)

    w = (E3|e2)
    w = w.expand()
    print('E3|e2 =',w)

    w = (E1|e1)
    w = (w.expand()).scalar()
    Esq = expand(Esq)
    print('%(E1\\cdot e1)/E^{2} =',simplify(w/Esq))

    w = (E2|e2)
    w = (w.expand()).scalar()
    print('%(E2\\cdot e2)/E^{2} =',simplify(w/Esq))

    w = (E3|e3)
    w = (w.expand()).scalar()
    print('%(E3\\cdot e3)/E^{2} =',simplify(w/Esq))
    return

# flip the current direction to follow convention
                    # that positive current flows from N1 to N2.
                    from lcapy import s

                    newelt = Element(I(elt.cpt.I), elt.nodes[1], elt.nodes[0],
                                     elt.name)
                    self.current_sources[key] = newelt
            else:
                raise ValueError('Unhandled element %s' % elt.name)


        A = self._A_matrix()
        Z = self._Z_vector()

        # Solve for the nodal voltages
        results = sym.simplify(A.inv() * Z);        

        # Create dictionary of node voltages
        self._V = Mdict({'0': Vs(0)})
        for n, node in enumerate(self.node_list[1:]):        
            self._V[node] = Vs(results[n])

        # Create dictionary of currents through elements
        self._I = {}
        for m, elt in enumerate(self.voltage_sources.values()):
            self._I[elt.name] = Is(results[m + self.num_nodes])

        for m, elt in enumerate(self.current_sources.values()):
            self._I[elt.name] = elt.cpt.I

        # Don't worry about currents due to initial conditions; these
        # are overwritten below.