derivative of MPolynomial_polydict

edit

In a for loop I generate a list of functions, then I take their derivatives with respect to different variables. For example, in the follwoing, Y[0,0] is a 'MPolynomial_polydict', and is equal to

x_0_1*x_0_1 + x_0_0*x_1_1 + x_0_2*x_0_2 + x_1_2*x_1_2 + x_0_0*x_2_2 + x_1_1*x_2_2 - 1.00000000000000

and list_of_variables[0] = x_0_1. When I take the derivative of the mentioned polynomial with respect to the mentioned variable I should get 2*x_0_1, but instead I get:

sage: derivative(Y[0,0],list_of_variables[0])
x_0_1

When I try this manually, it gives me the correct answer though:

sage: derivative(x_0_1*x_0_1 + x_0_0*x_1_1 + x_0_2*x_0_2 + x_1_2*x_1_2 + x_0_0*x_2_2 + x_1_1*x_2_2 - 1.00000000000000,x_0_1)
2*x_0_1

So, for some reason it considers the two x_0_1's in the first term of my polynomial as independent variables when I construct them automatically. Any ideas on why this problem arises, and how to fix this?

Edit:

This is how I build my vaiables:

names = [ [[] for i in range(n)] for j in range(n) ]

for i in range(n):
    for j in range(i+1):
        names[i][j] = (SR('x_' + str(j) + '_' + str(i)))
        names[j][i] = (SR('x_' + str(j) + '_' + str(i)))

R = PolynomialRing(CC,names[i][j] for i in range(n) for j in range(n)])
R.gens()
R.inject_variables()

More edit: Here is a the complete code:

#########################################################################
# Build variables, and the matrix corresponding to it
#########################################################################

def build_variables(n):
    names = [ [[] for i in range(n)] for j in range(n) ]

    for i in range(n):
        for j in range(i+1):
            names[i][j] = (SR('x_' + str(j) + '_' + str(i)))
            names[j][i] = (SR('x_' + str(j) + '_' + str(i)))

    return(names)


#########################################################################
# Define the function f that maps a matrix to the coefficients of its characteristic polynomial
#########################################################################
def CharPoly(Mat):
    X = matrix(Mat)
    n = X.ncols()
    C_X = X.characteristic_polynomial()
    Y = []
    for i in range(n):
        Y.append(C_X[i])
    return(Y)


############################################################################
# This solves that lambda problem
############################################################################
def lambda_siep(G,L,iter=100,epsilon = .1):
                        # G is any graph on n vertices
                        # L is the list of n desired distinct eigenvalues
                        # m is the number of itterations of the Newton's method
                        # epsilon: the off-diagonal entries will be equal to epsilon

    n = G.order()
    my_variables = build_variables(n)

    R = PolynomialRing(CC,[my_variables[i][j] for i in range(n) for j in range(n)])
    R.gens()
    R.inject_variables()


    X = [ [ R.gens()[n*i+j] for i in range(n) ] for j in range(n) ]

    Y = matrix(CharPoly(X)) - matrix(CharPoly(diagonal_matrix(L)))

    J = matrix(R,n)
    for i in range(n):
        for j in range(n):
            J[i,j] = derivative(Y[0][i],my_variables[j][j])

    B = diagonal_matrix(L) + epsilon * G.adjacency_matrix()

    count = 0
    while count < iter:

        T = [ B[i,j] for i in range(n) for j in range(n)]

        C = (J(T)).solve_right(Y(T).transpose())
        LC = list(C)

        B = B - diagonal_matrix([LC[i][0] for i in range(n)])

        count = count + 1

    return(B)