You should create the list res_coef from symbolic expressions, not directly from the coordinate functions representing the coefficients. To get the symbolic expression of the coordinate function res[ind], use the method expr(). Hence if you replace res[ind]**2 in your input In[5] by res[ind].expr()**2, everything works, as you can see on https://cloud.sagemath.com/projects/551a1e1d-9360-47bf-89ba-91603e96c7fe/files/Ask+sage.html