Bad matrix multiplication that gives the right answer
Python, N = 7
Using variable neighborhood search over elementary row and column operations.  Not every seed works quickly (or at all?), but seed 9 finds a solution for N = 7 in a few seconds.
$ time python3 search.py 7 9
[[0 3 6 6 1 6 6]
 [0 1 5 5 0 5 5]
 [1 3 1 2 1 8 1]
 [0 0 1 0 0 1 0]
 [1 7 9 9 0 9 9]
 [1 3 1 1 1 0 1]
 [0 0 1 1 0 1 1]]
[[5 5 1 1 0 4 0]
 [1 7 0 0 1 0 0]
 [0 0 9 4 0 1 5]
 [0 0 0 5 0 0 5]
 [2 2 1 1 7 4 0]
 [0 0 1 1 0 9 0]
 [0 2 0 0 0 0 0]]

real    0m4.909s
user    0m4.883s
sys     0m0.026s

Code
import itertools
import sys

import numpy as np

def search(depth, C, D, best):
    A = C + I
    B = D + 10 * I
    x = np.r_[A, B].flatten()
    score = (
        np.maximum(0, np.maximum(x - 9, 1 - x)).sum()
        - min((A == 0).sum(), -(-n * n // 2))
        - min((B == 0).sum(), -(-n * n // 2))
        + 2 * n
        - np.linalg.matrix_rank(A)
        - np.linalg.matrix_rank(B)
    )
    if best is None or score &lt; best:
        return C, D, score
    if depth == 0:
        return None

    tried = {}
    for i in reversed(range(n * (n - 1) * 5)):
        j = rng.integers(i + 1)
        mutation, tried[j] = tried.get(j, j), tried.get(i, i)
        p, mutation = divmod(mutation, (n - 1) * 5)
        q, mutation = divmod(mutation, 5)
        q = (p + q + 1) % n
        C1 = C.copy()
        D1 = D.copy()

        if mutation == 0:
            C1[p, :] += C1[q, :]
            D1[:, q] -= D1[:, p]
        elif mutation == 1:
            C1[p, :] -= C1[q, :]
            D1[:, q] += D1[:, p]
        elif mutation == 2:
            C1[:, p] += C1[:, q]
            D1[q, :] -= D1[p, :]
        elif mutation == 3:
            C1[:, p] -= C1[:, q]
            D1[q, :] += D1[p, :]
        elif p &lt; q:
            C1[[p, q], :] = C1[[q, p], :]
            D1[:, [p, q]] = D1[:, [q, p]]
        else:
            C1[:, [p, q]] = C1[:, [q, p]]
            D1[[p, q], :] = D1[[q, p], :]

        ret = search(depth - 1, C1, D1, best)
        if ret is not None:
            return ret

    return None

n, seed = map(int, sys.argv[1:])
rng = np.random.default_rng(seed)
I = np.eye(n, dtype=int)
score = None
k = rng.integers(0, 8, n)
C = np.diag(np.array([-10, -5, -2, -1, 1, 2, 5, 10])[k])
D = np.diag(np.array([-1, -2, -5, -10, 10, 5, 2, 1])[k])
depth = 1

while score != 0:
    ret = search(depth, C, D, score)
    if ret is not None:
        C, D, score = ret
        depth = 0
    depth += 1

A = C + I
B = D + 10 * I
assert (A @ B == 10 * A + B).all()
print(A)
print(B)