Matlab - Solving a system of 2 matricial equations : building each element of matrix solution

With Matlab, I have a matrix solution to find from 2 matricial equations (size matrix is 7x7). Here the 2 equations to solve with "`a`" and "`b`" are the unknow matrices and where `F1, F2, P1 and P2, D, D2, D` are known. Solving "a" and "b" would allow me to build a new matrix `P = a . P1 + b . P2` . (remark : `D` matrix is equal to : `D = a.a.D1 + b.b.D2` with D1and D2 diagonal matrices) :

1) a.a + a.P1.b.P2^T + b.P2.a.P1^T + b.b - Id = 0 (equation 1)

2) F1.a.P1 + F1.b.P2 + F2.a.P1 + F2.b.P2 − (a.P1 + b.P2).D = 0 (equation 2)

`Question 1)` Is there a solver in Matlab that would allow to directly find matrices "`a`" and "`b`" ?

`Question 2)` If there is no directly solving, I think that I have to solve each solution `a(i,j)` and `b(i,j`) by using a double loop on i and j indices and call fsolve in inner loop.

The problem is that I have to respect the order in matricial product and it is difficult to not make errors when I implement the code.

Here the 2 equations that I must solve on "`a`" and "`b`" unknown :

equation(1) :

equation 2):

So I tried to implement the solving of each matrix element, i.e a(i,j) and b(i,j)

For this, I am using a symbolic array `x[7 7 2]` and did :

    % eigenv_sp = P1
    % eigenv_xc = P2
    % FISH_sp = F1
    % FISH_xc = F2
    % FISH_eigen_sp = D1
    % FISH_eigen_xc = D2
    
    % Fisher matrices
    F1 = FISH_sp;
    F2 = FISH_xc;
    
    % Transposed matrix
    transp_eigenv_sp = eigenv_sp'
    transp_eigenv_xc = eigenv_xc'
    
    % Symbolic variables
    x = sym([7 7 2]);
    aP1 = sym([7 7]);
    aP1T = sym([7 7]);
    bP2 = sym([7 7]);
    bP2T = sym([7 7]);
    d1 = sym([7 7]);
    d2 = sym([7 7]);
    d = sym([7 7]);
    ad1 = sym([7 7]);
    bd2 = sym([7 7]);
    
    for k=1:7
      for l=1:7
        aP1(k,l) = sum(x(k,1:7,1).*eigenv_sp(1:7,l));
        bP2(k,l) = sum(x(k,1:7,2).*eigenv_xc(1:7,l));
    
        aP1T(k,l) = sum(x(k,1:7,1).*transp_egeinv_sp(1:7,l));
        bP2T(k,l) = sum(x(k,1:7,2).*transp_egeinv_xc(1:7,l));
    
        a2(k,l) = sum(x(k,1:7,1).*x(1:7,l,1));
        b2(k,l) = sum(x(k,1:7,2).*x(1:7,l,2));
    
        d1(k,l) = FkSH_ekgen_sp(k,l);
        d2(k,l) = FkSH_ekgen_xc(k,l);
        d(k,l) = d1(k,l) + d2(k,l);
    
        ad1(k,l) = sum(x(k,1:7,1).d1(1:7,l));
        bd2(k,l) = sum(x(k,1:7,2).d2(1:7,l));
      end
    end
    
    % Function that represents `(a,b)(i,j)` unknown element represented by `(y(i,j,1),y(i,j,2))`
    myfun=@(x,i,j) [
    
    % First equation
    sum(x(i,1:7,1).*x(1:7,j)) + sum(aP1(i,1:7).*bP2T(1:7,j)) + sum(bP2(i,1:7).*aP1T(1:7,j)) + sum(x(i,1:7).*x(1:7,2)) - eq(i,j);
    
    % second equation
    sum(F1(i,1:7).*aP1(1:7,j)) + sum(F1(i,1:7).*bP2(1:7,j)) + sum(F2(i,1:7).*aP1(1:7,j)) + sum(F2(i,1:7).*bP2(1:7,j)) - ...
    sum(aP1(i,1:7).*ad1(1:7,j)) + sum(bP2(i,1:7).*ad1(1:7,j)) + sum(adP1(i,1:7).*dP2(1:7,j)) + sum(bP2(i,1:7).*dP2(1:7,j));
    ]
    
    % Solution of system of non linear equations with loop
    y = zeros(7, 7, 2);
    for i=1:7
      for j=1:7
        a0 = 1e7;
        b0 = 1e7;
        x0 = [ a0 b0];
        y(i,j,:) = fsolve(@(x)myfun(x,i,j),x0)
        end
    end

But at the execution, I have the following error :

    Index exceeds matrix dimensions.
    
    Error in sym/subsref (line 814)
                R_tilde = builtin('subsref',L_tilde,Idx);
    
    Error in compute_solving_Matricial_Equations_dev (line 60)
        aP1(k,l) = sum(x(k,1:7,1).*eigenv_sp(1:7,l));

I don't understand where is the issue, I only iterate on index `1:7`, if someone could see what's wrong ?

UPDATE 1: I wonder if I have the right to use for a given index `"i" the symbolic expression `x(1:7,i,1)` whereas this array is not known numerically but defined only from a symbolic point of view : do you understand what I mean ?

UPDATE 2 : Maybe I have a clue.

If I do an expand of the symbolic array `x` (of size `(7,7,2)`), I get :

    >> expand(sum(x(1,1:7).*x(1:7,1)))
    Index exceeds matrix dimensions.
    
    Error in sym/subsref (line 814)
                R_tilde = builtin('subsref',L_tilde,Idx);

So, the only thing to done I think is to explicitly write the 7x7 equations

with 7x7 = 49 couples unknown `x(i,j,1)` and `x(i,j,2)` which corresponds to `a(i,j)` and `b(i,j)` matrices.

UPDATE 3: I think the key point is to put into function `myfun` the 2 matrical equations but with taking into account of the 49 (7x7) couples `(a(i,j) and b(i,j))` unknown that I defined previously as symbolic variables with array `x[7 7 2]`. And Second point, I must to "make understand" to Matlab with `fsolve` that we are searching for each iteration (double loop on `(i,j)`) the value of `a(i,j)` and `b(i,j)` given the whole 49 equations. Maybe it would be necessary to expand for example a symbolic expression of `x(i,1:7)` or `x(1:7,j)` in `myfun` function but I don't know how to proceed.

Any help is welcome, thanks in advance.