Claro que se puede, con Python el cielo es el límite... XD.

Bromas aparte, no especificas que suma de Riemann quieres aplicar (derecha, izquierda, máximo, mínimo). Voy a usar la suma izquierda de acuerdo a tu ejemplo:

Para la gráfica uso <code>[matplotlib][1]</code> y <code>[NumPy][2]</code> para los arrays:

    import numpy as np
    import matplotlib.pyplot as plt


    def riemannplot(f, a, b, ra, rb, n):
        # f es la función 
        # a y b son los limites del eje x para graficar la funcion f
        # ra y rb son los limites del intervalo en el eje x del que queremos calcular la suma
        # n es el numero de rectangulos que calcularemos
    
        atenuacion = (b-a)/100
        x = np.arange(a, b+atenuacion, atenuacion)
        plt.plot(x, f(x), color='green')
    
        delta_x = (rb-ra)/n
        riemannx = np.arange(ra, rb, delta_x)
        riemanny = f(riemannx)
        riemann_sum = sum(riemanny*delta_x)

        plt.bar(riemannx,riemanny,width=delta_x,alpha=0.5,facecolor='orange')
        plt.xlabel('x')
        plt.ylabel('f(x)')
        plt.title('Suma izquierda de Riemann para f(x)')
        area_text = plt.figtext(0.1,0.02, "Suma de Riemann: " + str(riemann_sum), color='r')
        plt.show()



    def f(x):
        return x**2

    riemannplot(f, 0, 1.1, 0, 1, 10)


No se si tienes que hacerlo mediante <code>[SymPy][3]</code> obligatoriamente o no.
La gráfica que nos crea es esta:

[![introducir la descripción de la imagen aquí][4]][4]


  [1]: http://matplotlib.org/
  [2]: http://www.numpy.org/
  [3]: http://www.sympy.org/es/
  [4]: https://i.sstatic.net/9OOTN.png