Integrate Legendre polynomials in SymPy and use these integrals as coefficients

I am trying to make a simple example in SymPy to compute some coefficients and then use them in a sum of legendre polynomials. Finally plot it. Very simple but can not make it work. I want to use it for my electromagnetism course. I get errors in both the attempts below:

%matplotlib inline
from sympy import *
x, y, z = symbols('x y z')
k, m, n = symbols('k m n', integer=True)
f, step, potential = symbols('f step potential', cls=Function)
var('n x')
A=SeqFormula(2*(2*m+1)*Integral(legendre(2*m+1,x),(x,0,1)).doit(),(m,0,oo)).doit()
Sum(A.coeff(m).doit()*legendre(2*m+1,x),(m,0,10)).doit()
B=Tuple.fromiter(2*(2*n+1)*Integral(legendre(2*n+1,x),(x,0,1)).doit() for n in range(50))
Sum(B[m]*legendre(2*m+1,x),(m,0,10)).doit()

Here is a part of an script in Mathematica of what I would like to replicate:

Nn = 50; 
Array[A, Nn] 
For[i = 0, i <= Nn, i++, A[i + 1] = Integrate[LegendreP[2*i + 1, x]*(2*(2*i + 1) + 1), {x, 0, 1}]]; 
Step = Sum[A[n + 1]*LegendreP[2*n + 1, #], {n, 0, Nn}] & Plot[Step[x], {x, -1, 1}]

Copyright Notice：Content Author:「Oscar Reula」，Reproduced under the CC 4.0 BY-SA copyright license with a link to the original source and this disclaimer.
Link to original article：https://stackoverflow.com/questions/49259905/integrate-legendre-polynomials-in-sympy-and-use-these-integrals-as-coefficients