You might get close enough with jacobian to get the full output:
syms x y real
% simple 2 x 2 matrix:
M = [sin(x+y) cos(x*y);-cos(x/(x^2+y^2)) tan(x)];
% Should return this:
% M =
%
% [ sin(x + y), cos(x*y)]
% [-cos(x/(x^2 + y^2)), tan(x)]
%
% This will turn the 2x2 matrix M into a 4x1 column-array, which jacobian
% can handle:
J = jacobian(M(:),[x,y])
% should return
% J =
%
% [ cos(x + y), cos(x + y)]
% [sin(x/(x^2 + y^2))*(1/(x^2 + y^2) - (2*x^2)/(x^2 + y^2)^2), -(2*x*y*sin(x/(x^2 + y^2)))/(x^2 + y^2)^2]
% [ -y*sin(x*y), -x*sin(x*y)]
% [ tan(x)^2 + 1, 0]
Here you get all 4 x 2 derivatives, you will have to re-sort the rows into the ordering you expect. For this toy-model-example I get:
J3D = reshape(J,[2,2,2]);
to produce something like what you describe. In your case this might be a rather big array of the size n1 x n2 x n_vars - if I understood your problem right.
HTH
