def element_matrix(phi, Omega_e, symbolic=True):
n = len(phi)
A_e = sp.zeros((n, n))
X = sp.Symbol('X')
if symbolic:
h = sp.Symbol('h')
else:
h = Omega_e[1] - Omega_e[0]
detJ = h/2 # dx/dX
for r in range(n):
for s in range(r, n):
A_e[r,s] = sp.integrate(phi[r]*phi[s]*detJ, (X, -1, 1))
A_e[s,r] = A_e[r,s]
return A_e