matrix(
    [W, rho * v_x * W^3, rho * v_y * W^3, rho * v_z * W^3, 0],
    [h * W^2 * v_x, rho * h * W^2 * (1 + 2 * W^2 * v_x^2), 2 * rho * h * W^4 * v_x * v_y, 2 * rho * h * W^4 * v_x * v_z, rho * W^2 * v_x * gamma],
    [h * W^2 * v_y, 2 * rho * h * W^4 * v_x * v_y, rho * h * W^2 * (1 + 2 * W^2 * v_y^2), 2 * rho * h * W^4 * v_y * v_z, rho * W^2 * v_y * gamma],
    [h * W^2 * v_z, 2 * rho * h * W^4 * v_x * v_z, 2 * rho * h * W^4 * v_y * v_z, rho * h * W^2 * (1 + 2 * W^2 * v_z^2), rho * W^2 * v_z * gamma],
    [-(W - epsilon + epsilon * gamma - h * W^2), -rho * W^3 * v_x * (1 - 2 * h * W), -rho * W^3 * v_y * (1 - 2 * h * W), -rho * W^3 * v_z * (1 - 2 * h * W), rho * (1 + W^2 * gamma - gamma)] 
);