Python et méthode ADI avec matrices

import numpy as np
import scipy as sp
import scipy.linalg as la

n=3
h = 1.0/(n+1)
A = np.zeros((n,))
A[:2] = (2.0/h**2, -1.0/h**2)

A = la.toeplitz(A)

AA = np.identity(n)

B = np.kron(AA,A)

I = np.identity(n)

II = np.zeros((n,))
II[:2] = (2.0/h**2, -1.0/h**2)
II = la.toeplitz(II)

C= np.kron(II,I)

J=np.identity(9)
def ADI(B,C,b,Imax,eps,r,x0,y0):
    x=x0
    y=y0
    i=0
    while (i<Imax):
        y=la.solve((np.dot(r,J)+b),np.dot((np.dot(r,J)-C),x)+b)
        x=la.solve((np.dot(r,J)+C),np.dot((np.dot(r,J)-B),y)+b)
        i=i+1
    return (x,y,i)

b=np.ones((9,1))
x0=np.ones((9,1))
y0=np.ones((9,1))
print(ADI(B,C,b,1000,10**(-3),1,x0,y0))