How to update the results from a solve function to the original L and U matrices?

I can find each equation to solve but I cannot seem to update each computation's result in the original matrix.

I tried to itemset in the matrix, replacing the u variables with actual result from the solve function. But then it gives the following error: TypeError: can't multiply sequence by non-int of type 'Symbol'

Dont know what to do else!

from sympy import *
import numpy as np


#input matrix

a = np.array([[1,3,5],\
              [2,4,7],\
              [1,1,0]])


d = len(a)
d1 = d-1
d2 = (d1*(d1+1))/2
d3 = (d*(d+1))/2


symbols_dict = dict(('L%d'%k, symbols('L%d'%k)) for k in range(d2))
locals().update(symbols_dict)
symbols_dict = dict(('U%d'%k, symbols('U%d'%k)) for k in range(d3))
locals().update(symbols_dict)


# L decomposition

b = np.array([[1,0,0],\
              [L0,1,0],\
              [L1,L2,1]])

# U decomposition

c = np.array([[U0,U1,U2],\
              [0,U3,U4],\
              [0,0,U5]])


t=[[[]for i in range(d)]for j in range(d)]
s=[[[]for i in range(d)]for j in range(d)]

for j in range(d):
    for i in range(d):
        m=0
        for k in range(d):
            m = m+(b[j,k]*c[k,i])
            t[j][i] = m
        s[j][i]= solve(t[j][i]-a[j,i])
        c.itemset((j,i),s[j][i])

It should give me the numerical results for each equation solved.

The corrected code is as follows:

L0,L1,L2 = symbols('L0 L1 L2')

U0,U1,U2,U3,U4,U5 = symbols('U0 U1 U2 U3 U4 U5')
.
.
.


for j in range(d):
    for i in range(d):
        m=0
        for k in range(d):
            m = m+(b[j,k]*c[k,i])
            s= solve(m-a[j,i])
            if j>i:
                b[j,i]=s
            else:
                c[j,i]=s