1.

using TyOptimization
fun = x-> begin
F1 = x[1]^2 - x[2] - 2
F2 = x[2]^2 - 2 * x[1] - 4
return [F1,F2]
end

x0 = [2, 2]
x, = fsolve(fun, x0);

2.

rho(z, y, x) = 2 + z
xmin = -1
xmax = 1
ymin(x) = -sqrt(1 - x^2)
ymax(x) = sqrt(1 - x^2)
zmin(x, y) = sqrt(x^2 + y^2)
zmax(x, y) = 4 - x^2 - y^2
m, e = integral3(rho, xmin, xmax, ymin, ymax, zmin, zmax)