In [2]:
import sympy as sp
sp.init_printing()
from IPython.display import display
In [61]:
x,y,z = sp.symbols('x, y, z')                # set symbols

ex = x**2 + y**2 + z                         # create an expression

display(ex)                                  # display expression (better than print)

ex2 = ex.subs(y+1,x)                         # substitute y+1 for x
ex3 = ex.subs( [(y+1,x), (z, 7)])            # multiple substitution

ex.evalf(subs={x:3, y:4, z:sp.pi})           # numerical evaluation
ex.evalf(100, subs={x:3, y:4, z:sp.pi})      # 100 digits of accuracy

ex = (x**2 + 3*x + 2)/(x+1)
ex.simplify()                                # simplify: three versions
sp.simplify(ex)
sp.simplify( (x**2 + 3*x + 2)/(x+1) )

ex = (x+2)*(x-3)                             # expand : three versions
ex.expand()
sp.expand(ex)
sp.expand( (x+2)*(x-3) )

ex = x**3 - x**2 + x - 1
ex.factor()
sp.factor(ex)
sp.factor( x**3 - x**2 + x -1 )

# collect, cancel, apart, trigsimp, expand_trig, powsimp
# expand_power_exp, expand_power_base, pow_denest, expand_log, etc.

ex = sp.exp(x*y*z)
ex.diff(x)                   # differentiate: three versions
sp.diff(ex,x)
sp.diff(sp.exp(x**2), x)
sp.diff(ex, x,x,x)           # 3rd derivative
ex.diff(x,y,z)               # 3rd mixed derivative

exD = sp.Derivative(ex, x,y,z)  # expression with derivative
exD.doit()                      # evaluate the derivative expression

ex = x**2 + y**2
ex.integrate(x)                      # integrate: three versions
sp.integrate(ex, x)                  # (note, the constant is left off)
sp.integrate(x**2, x)      
sp.integrate(ex, (x,0,sp.oo))        # definite integral
sp.integrate(ex, (x,1,2), (y, 0, z)) # double integral

exI = sp.Integral(ex, x)             # expression with integral
exI = sp.Integral(ex, (x,0,2))       # include bounds
exI.doit()                           # evaluate the integral expression

ex = sp.sin(x)/x         
exL = sp.limit(ex, x, 0)             # limit
exL = sp.Limit(ex, x, 0)             # expression with limit
exL.doit()                           # evaluate the limit expression

sp.solve(x**2-y, x)                  # solve
ex = sp.Eq(x**2, y)                  # create an equality: x**2 = y
sp.solve(ex, x)
sp.solve( (x-y+2, x+y-3), (x, y) )   # solve system of equations
sp.solve( (x*y-7, x+y-6), (x, y) )   # nonlinear: try replace 6 with z

f = sp.Function('f')                 # function variable: 2 ways
f = sp.symbols('f', cls=sp.Function)

diffEq = f(x).diff(x,x) - 2*f(x).diff(x) + f(x) - sp.sin(x)
sp.dsolve(diffEq, f(x))              # solve ODE for f(x)

m11, m12, m21, m22 = sp.symbols("m11, m12, m21, m22")
b1, b2             = sp.symbols("b1, b2")
A = sp.Matrix([ [m11, m12], [m21, m22] ])
b = sp.Matrix([b1,b2])
A.inv()*b      
A**-1 * b
$$x^{2} + y^{2} + z$$
Out[61]:
$$\left[\begin{matrix}b_{1} \left(\frac{1}{m_{11}} + \frac{m_{12} m_{21}}{m_{11}^{2} \left(m_{22} - \frac{m_{12} m_{21}}{m_{11}}\right)}\right) - \frac{b_{2} m_{12}}{m_{11} \left(m_{22} - \frac{m_{12} m_{21}}{m_{11}}\right)}\\- \frac{b_{1} m_{21}}{m_{11} \left(m_{22} - \frac{m_{12} m_{21}}{m_{11}}\right)} + \frac{b_{2}}{m_{22} - \frac{m_{12} m_{21}}{m_{11}}}\end{matrix}\right]$$