The goal of the matlab program is to find the following optimization problem:
The optimized variable are four symmetric matrices A0, A1, B0, B1 of dimension 4. rho is a given 16 * 16 matrix (assume it is identity matrix).
The constraint is -I <=A0, A1, B0, B1<= I, here A<=B means B-A is semidefinite positive.
The objective function is max Tr[(A0 otimes B0 + A1 otimes B0 + A0 otimes B1 – A1 otimes B1)*rho], where otimes means the Kronecker product.

I write the following matlab program:

A0 = sym('A0', [4 4]);
B0 = sym('B0', [4 4]);
A1 = sym('A1', [4 4]);
B1 = sym('B1', [4 4]);
rho= eye(16);
constraints = {
eye(4) + A0 >= 0,
eye(4) - A0 >= 0, 
eye(4) - B0 >= 0,
eye(4) - B0 >= 0,
eye(4) + A1 >= 0,
eye(4) - A1 >= 0, 
eye(4) - B1 >= 0,
eye(4) - B1 >= 0,
};
D00 = kron(A0, B0);
D01 = kron(A0, B1);
D10 = kron(A1, B0);
D11 = kron(A1, B1);
objective = trace(rho * (D00 + D01 + D10 - D11));

However I don’t know how to do the final solution step because here the variables are matrices so reshape is needed. With the hint of chatgpt, it says

options = optimoptions('fmincon', 'Display', 'iter', 'MaxFunctionEvaluations', 10000);
[x, fval] = fmincon(@(x) -double(subs(objective, {A0, B0, A1, B1}, {reshape(x(1:16), 4, 4), reshape(x(17:32), 4, 4), reshape(x(33:48), 4, 4), reshape(x(49:64), 4, 4)})), ...
rand(64, 1), [], [], [], [], [], [], @(x)deal(constraints), options);

but it reports error: Entries in second argument must be scalar.