Solving system of 9 nonlinear equaitons in 16 variables

I have a system of equations as follows:

I am not able to use fsolve as it says in the documentaiton that the number of variables should be as same as the number of equations. I found this on the MathWorks which says that it can be done with fsolve. Please let me know if it can be solved by any other method or by using fsolve. It will also suffice if I can know the solution exists.

I am writing the MATLAB code that I have written using fsolve.

f = @(x) [x(1)*x(9)  + x(2)*x(12) + x(3)*x(15)  - 13;
          x(1)*x(10) + x(2)*x(13) + x(3)*x(16)  - 15;
          x(1)*x(11) + x(2)*x(14) - x(3)*(x(9)  + x(13)) + 1;
          x(4)*x(9)  + x(5)*x(12) + x(6)*x(15)  - 9;
          x(4)*x(10) + x(5)*x(13) + x(6)*x(16)  - 24;
          x(4)*x(11) + x(5)*x(14) - x(6)*(x(9)  + x(13));
          x(7)*x(9)  + x(8)*x(12) - x(15)*(x(1) + x(5)) - 7;
          x(7)*x(10) + x(8)*x(13) - x(16)*(x(1) + x(5)) -2;
          x(7)*x(11) + x(8)*x(14) + (x(1)+x(5))*(x(9)+x(13)) - 35];

A = zeros(1,9);
fsolve(f, A)

Answer 1

[Torsten]

x0 = -10*ones(16,1);
AB = [13 15 -1;9 24 0;7 2 35];
options = optimset('TolFun',1e-16,'TolX',1e-16);
x = fmincon(@(x)fun(x,AB),x0,[],[],[],[],[],[],[],options);

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

A = [x(1) x(2) x(3);x(4) x(5) x(6);x(7) x(8) -(x(1)+x(5))]

A = 3×3

-0.2352 1.4168 -3.4167 -4.3530 2.0720 -3.1488 -0.9938 -5.0327 -1.8368

B = [x(9) x(10) x(11);x(12) x(13) x(14);x(15) x(16) -(x(9)+x(13))]

B = 3×3

0.7008 -2.2052 -1.3224 -0.1068 1.3766 -5.9351 -3.8974 -3.6675 -2.0774

A*B-AB

ans = 3×3

1.0e-06 * -0.0907 -0.1745 -0.0038 -0.0097 -0.1094 0.0113 -0.0376 0.0691 -0.1961

function obj = fun(x,AB)
  A = [x(1) x(2) x(3);x(4) x(5) x(6);x(7) x(8) -(x(1)+x(5))];
  B = [x(9) x(10) x(11);x(12) x(13) x(14);x(15) x(16) -(x(9)+x(13))];
  M = A*B - AB;
  M = M(:);
  obj = sum(M.^2);
end