One approach:
x = [1,2,3,4,5,6,7];
y = [1,2,3,4,5,6,7];
[X,Y] = meshgrid(x,y);
xv = X(:);
yv = Y(:);
% -3*w*sin(x)+ 5*q*sin(y)-10 = 0
% 3*w*cos(x) - 5*q*cos(y) = 0
FM = @(w,q,x,y) [-3*w*sin(x) 5*q*sin(y)-10; 3*w*cos(x) -5*q*cos(y)];
for k = 1:numel(xv)
B(:,k) = fsolve(@(b) FM(b(1),b(2),xv(k),yv(k)), [100;100]);
end
Here ‘B(1,:)’ is ‘w’, and ‘B(2,:)’ is ‘q’. They are due to the matching ‘xv(k)’ ‘yv(k)’ combineations for each ‘k’.
Experiment to get the result you want.
