def MSI(A,b,eps):
    if not np.all(np.linalg.eigvals(A) > 0):
        T = A
        A = np.transpose(T).dot(A)
        b = np.transpose(T).dot(A)

    mu = 2/(np.linalg.norm(A) + eps)
    B = np.eye(A.shape[0]) - mu * A
    c = mu*b
    xcur = np.random.random(A.shape[0])
    xcur = np.reshape(xcur, (3,1))
    xnext = B@xcur + c


    while (np.linalg.norm(xcur - xnext) > eps):
        xcur = xnext
        xnext = B.dot(xcur) + c
    return xnext