こういう事ではないかと思います。
G = [10 11 ; 12 13];
H = [-2 -1 ; -2 -3];
% f = @(x)[x(1) .*x(2) -5*x(2) +G(:,:), x(1) .^3-x(2) .^2+H(:,:)];
f = @(x)[x(:,1:2).*x(:,3:4)-5*x(:,3:4)+G(:,:), x(:,1:2).^3-x(:,3:4).^2+H(:,:)];
I = [1 1 ; 1 1];
J = [1 1 ; 1 1];
x0 = [I J]; % 初期値
x = fsolve(f,x0) % ,options)