syms x y z;

f(x,y,z) = [x*y^3; x*y*z; y^2*z^4];

disp(jacobian(f));