How to use the sympy.symbols function in sympy

print('B =', B2d)
    print('A + B =', A2d + B2d)
    print('A - B =', A2d - B2d)

    # TODO: add this back when we drop Sympy 1.3. The 64kB of output is far too
    # printer-dependent
    if False:
        print('AB =', A2d * B2d)

    a = g2d.mv('a','vector')
    b = g2d.mv('b','vector')

    print(r'a|\f{\overline{A}}{b}-b|\f{\underline{A}}{a} =',((a|A2d.adj()(b))-(b|A2d(a))).simplify())

    m4d = Ga('e_t e_x e_y e_z', g=[1, -1, -1, -1],coords=symbols('t,x,y,z',real=True))

    T = m4d.lt('T')

    print('g =', m4d.g)

    print(r'\underline{T} =',T)
    print(r'\overline{T} =',T.adj())

    print(r'\f{\det}{\underline{T}} =',T.det())
    print(r'\f{\mbox{tr}}{\underline{T}} =',T.tr())

    a = m4d.mv('a','vector')
    b = m4d.mv('b','vector')

    print(r'a|\f{\overline{T}}{b}-b|\f{\underline{T}}{a} =',((a|T.adj()(b))-(b|T(a))).simplify())

def Distorted_manifold_with_scalar_function():
    coords = symbols('x y z')
    (ex,ey,ez,grad) = MV.setup('ex ey ez',metric='[1,1,1]',coords=coords)
    mfvar = (u,v) = symbols('u v')
    X = 2*u*ex+2*v*ey+(u**3+v**3/2)*ez
    MF = Manifold(X,mfvar,I=MV.I)

    (eu,ev) = MF.Basis()

    g = (v+1)*log(u)
    dg = MF.Grad(g)
    print 'g =',g
    print 'dg =',dg
    print 'dg(1,0) =',dg.subs({u:1,v:0})
    G = u*eu+v*ev
    dG = MF.Grad(G)
    print 'G =',G
    print 'P(G) =',MF.Proj(G)
    print 'zcoef =',simplify(2*(u**2 + v**2)*(-4*u**2 - 4*v**2 - 1))
    print 'dG =',dG

F = E+I*B

    print('B = \\bm{B\\gamma_{t}} =',B)
    print('E = \\bm{E\\gamma_{t}} =',E)
    print('F = E+IB =',F)
    print('J =',J)
    gradF = m4d.grad*F
    gradF.Fmt(3,'grad*F')

    print('grad*F = J')
    (gradF.get_grade(1)-J).Fmt(3,'%\\grade{\\nabla F}_{1} -J = 0')
    (gradF.get_grade(3)).Fmt(3,'%\\grade{\\nabla F}_{3} = 0')

    (alpha,beta,gamma) = symbols('alpha beta gamma')

    (x,t,xp,tp) = symbols("x t x' t'")
    m2d = Ga('gamma*t|x',g=[1,-1])
    (g0,g1) = m2d.mv()

    R = cosh(alpha/2)+sinh(alpha/2)*(g0^g1)
    X = t*g0+x*g1
    Xp = tp*g0+xp*g1
    print('R =',R)

    print(r"#%t\bm{\gamma_{t}}+x\bm{\gamma_{x}} = t'\bm{\gamma'_{t}}+x'\bm{\gamma'_{x}} = R\lp t'\bm{\gamma_{t}}+x'\bm{\gamma_{x}}\rp R^{\dagger}")

    Xpp = R*Xp*R.rev()
    Xpp = Xpp.collect()
    Xpp = Xpp.trigsimp()
    print(r"%t\bm{\gamma_{t}}+x\bm{\gamma_{x}} =",Xpp)
    Xpp = Xpp.subs({sinh(alpha):gamma*beta,cosh(alpha):gamma})

def get_CC_operators():
    """
    Returns a tuple (T1,T2) of unique operators.
    """
    i = symbols('i', below_fermi=True, cls=Dummy)
    a = symbols('a', above_fermi=True, cls=Dummy)
    t_ai = AntiSymmetricTensor('t', (a,), (i,))
    ai = NO(Fd(a)*F(i))
    i, j = symbols('i,j', below_fermi=True, cls=Dummy)
    a, b = symbols('a,b', above_fermi=True, cls=Dummy)
    t_abij = AntiSymmetricTensor('t', (a, b), (i, j))
    abji = NO(Fd(a)*Fd(b)*F(j)*F(i))

    T1 = t_ai*ai
    T2 = Rational(1, 4)*t_abij*abji
    return (T1, T2)

def expr_body(expr, **kwargs):
    if not hasattr(expr, '__len__'):
        # Defined in terms of some coordinates
        xyz = set(sp.symbols('x[0], x[1], x[2]'))
        xyz_used = xyz & expr.free_symbols
        assert xyz_used <= xyz
        # Expression params which need default values
        params = (expr.free_symbols - xyz_used) & set(kwargs.keys())
        # Body
        expr = ccode(expr).replace('M_PI', 'pi')
        # Default to zero
        kwargs.update(dict((str(p), 0.) for p in params))
        # Convert
        return expr
    # Vectors, Matrices as iterables of expressions
    else:
        return [expr_body(e, **kwargs) for e in expr]

def MV_setup_options():
    (e1,e2,e3) = MV.setup('e_1 e_2 e_3','[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    (e1,e2,e3) = MV.setup('e*1|2|3','[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    (e1,e2,e3) = MV.setup('e*x|y|z','[1,1,1]')
    v = MV('v', 'vector')
    print(v)

    coords = symbols('x y z')
    (e1,e2,e3,grad) = MV.setup('e','[1,1,1]',coords=coords)
    v = MV('v', 'vector')
    print(v)

    return

# (that means the dimensions don't match)
                j += 1
                
        # set system dimensions
        self.n = n
        self.m = m

        logging.debug("--> system: {}".format(n))
        logging.debug("--> input : {}".format(m))

        # next, we look for integrator chains
        logging.debug("Looking for integrator chains")

        # create symbolic variables to find integrator chains
        x_sym = ([sp.symbols('x%d' % k, type=float) for k in xrange(1,n+1)])
        u_sym = ([sp.symbols('u%d' % k, type=float) for k in xrange(1,m+1)])

        fi = self.ff_sym(x_sym, u_sym)

        chaindict = {}
        for i in xrange(len(fi)):
            # substitution because of sympy difference betw. 1.0 and 1
            if isinstance(fi[i], sp.Basic):
                fi[i] = fi[i].subs(1.0, 1)

            for xx in x_sym:
                if fi[i] == xx:
                    chaindict[xx] = x_sym[i]

            for uu in u_sym:
                if fi[i] == uu:
                    chaindict[uu] = x_sym[i]

# -*- coding: utf-8 -*-
"""Exercise 10.15 from Kane 1985."""

from __future__ import division
from sympy import expand, sin, cos, solve, symbols, trigsimp, together
from sympy.physics.mechanics import ReferenceFrame, Point, Particle
from sympy.physics.mechanics import dot, dynamicsymbols, msprint
from util import generalized_inertia_forces, generalized_inertia_forces_K
from util import partial_velocities


q1, q2, q3 = q = dynamicsymbols('q1:4')
q1d, q2d, q3d = qd = dynamicsymbols('q1:4', level=1)
u1, u2, u3 = u = dynamicsymbols('u1:4')
g, m, L, t = symbols('g m L t')
Q, R, S = symbols('Q R S')

# reference frame
N = ReferenceFrame('N')

# points and velocities
# Simplify the system to 7 points, where each point is the aggregations of
# rods that are parallel horizontally.
pO = Point('O')
pO.set_vel(N, 0)

pP1 = pO.locatenew('P1', L/2*(cos(q1)*N.x + sin(q1)*N.y))
pP2 = pP1.locatenew('P2', L/2*(cos(q1)*N.x + sin(q1)*N.y))
pP3 = pP2.locatenew('P3', L/2*(cos(q2)*N.x - sin(q2)*N.y))
pP4 = pP3.locatenew('P4', L/2*(cos(q2)*N.x - sin(q2)*N.y))
pP5 = pP4.locatenew('P5', L/2*(cos(q3)*N.x + sin(q3)*N.y))
pP6 = pP5.locatenew('P6', L/2*(cos(q3)*N.x + sin(q3)*N.y))

import importlib
from sympy import symbols, cos, sin, chebyshevt
import numpy as np
from mpi4py import MPI
from shenfun import inner, div, grad, TestFunction, TrialFunction, Array, \
    Function, TensorProductSpace, Basis, extract_bc_matrices

comm = MPI.COMM_WORLD

# Collect basis and solver from either Chebyshev or Legendre submodules
family = sys.argv[-1].lower() if len(sys.argv) == 2 else 'chebyshev'
base = importlib.import_module('.'.join(('shenfun', family)))
BiharmonicSolver = base.la.Biharmonic

# Use sympy to compute a rhs, given an analytical solution
x, y = symbols("x,y")
a = 1
b = -1
if family == 'jacobi':
    a = 0
    b = 0
ue = (sin(2*np.pi*x)*cos(2*y))*(1-x**2) + a*(0.5-9./16.*x+1./16.*chebyshevt(3, x)) + b*(0.5+9./16.*x-1./16.*chebyshevt(3, x))
#ue = (sin(2*np.pi*x)*cos(2*y))*(1-x**2) + a*(0.5-0.6*x+1/10*legendre(3, x)) + b*(0.5+0.6*x-1./10.*legendre(3, x))
fe = ue.diff(x, 4) + ue.diff(y, 4) + 2*ue.diff(x, 2, y, 2)

# Size of discretization
N = (30, 30)

if family == 'chebyshev':
    assert N[0] % 2 == 0, "Biharmonic solver only implemented for even numbers"

#SD = Basis(N[0], family=family, bc='Biharmonic')

def derivatives_in_spherical_coordinates():
    Print_Function()
    X = (r,th,phi) = symbols('r theta phi')
    curv = [[r*cos(phi)*sin(th),r*sin(phi)*sin(th),r*cos(th)],[1,r,r*sin(th)]]
    (er,eth,ephi,grad) = MV.setup('e_r e_theta e_phi',metric='[1,1,1]',coords=X,curv=curv)

    f = MV('f','scalar',fct=True)
    A = MV('A','vector',fct=True)
    B = MV('B','grade2',fct=True)

    print('f =',f)
    print('A =',A)
    print('B =',B)

    print('grad*f =',grad*f)
    print('grad|A =',grad|A)
    print('-I*(grad^A) =',-MV.I*(grad^A))
    print('grad^B =',grad^B)
    return