Suppose $X_1 \sim f_{X_1}(x_1)$, $X_2 \sim f_{X_2}(x_2)$ are random variables with known probability density function.
Is there any way to compute the probability density function of a bivariate function $g(x_1,x_2)$, assuming on specific finite support $x_1 \in [a_1,b_1]$ and $x_2 \in [a_2,b_2]$ without Monte Carlo sampling or any other sampling method? $X_1$ and $X_2$ are independent, and $g $ is "smooth".
I'm looking for a way to possibly implement the computation numerically for an arbitrary $g$.
HINT
Let's rename your variables $X,Y$ for easier notation. So you are defining $$ Z = g(X,Y) $$ and asking what is the probability density function of $Z$. For simplicity, let's assume a particular (very simple) $g$, let's say $Z = X+Y$.
$Z=z$ means $X+Y=z$, so $X=z-Y$; let's condition on $Y=y$, then you end up with $$ f_Z(z) = \int_\mathbb{R} f_X(z-y) f_Y(y) dy. $$
So if $g(x,y)$ has a nice form allowing to solve $z = g(x,y)$ for one of the variables in terms of the other, you can use this technique to get the pdf...