I am trying to derive the following equation from Wikipedia on probability density function which is that for functions that are not monotonic, the probability density for $y$ is:
$$f_Y(y) = \sum \limits_{k=1}^{n(y)} \left| \frac{\mathrm{d}}{\mathrm{d}y} g^{-1}|_k(y) \right| f_X (g^{-1}|_k(y)) .$$
My attempt is the following: first of all, I am not working with axioms but just with real analysis and some intuition about what probability should be. Let us assume that we know $f_X(x)$ and we want to find out $f_Y(y)$. Here, $x \in [x_1, x_2)$ and $y \in [y_1, y_2)$. Also, assume that the map between variables is $g: x \mapsto y$. Now, assume that $[x_1, x_2)$ can be divided into finite disjoint intervals such that their union is $[x_1, x_2)$. Call these intervals $U_k = [x_k, x_{k+1})$ for $k=1,2, \ldots, N$. Assume that in each interval $g$ restricted to $U_k$ (call it $g|_k$) is a bijection.
\begin{align} \mathrm{Prob}_Y ([a,b]) & = \operatorname{Prob}_X(g^{-1}([a,b])) \\[8pt] & = \operatorname{Prob}_X \left(\bigcup_k (g^{-1}([a,b]) \cap U_k)\right) \\[8pt] & = \sum_k \operatorname{Prob}_X (g^{-1}([a,b]) \cap U_k) \end{align}
Now, because I model probability using continuous distributions, I can translate this result by using change of variables in each segment as there function is bijection.
$$ \operatorname{Prob}_X (g^{-1}([a,b]) \cap U_k) = \int_{g^{-1}([a,b]) \cap U_k} f_X(x) \, \mathrm{d}x = \int_{g|_k (g^{-1}([a,b]) \cap U_k)} f_X(g^{-1}|_k(y)) \left| g^{-1}|_k' \right|(y) \, \mathrm{d}y$$
Now, we combine all computations to have the following result.
$$ \int \limits_a^b f_Y(y) \, \mathrm{d}y = \sum_k \int_{g|_k (g^{-1}([a,b]) \cap U_k)} f_X(g^{-1}|_k(y)) \left| g^{-1}|_k' \right|(y) \, \mathrm{d}y$$
I don't understand how to make the next step. I guess I would have to differentiate both sides with respect to $b$ and use fundamental theorem of calculus on the left hand side, but I am not sure how to differentiate (if possible) the right side.
I think you've basically got it. Perhaps writing it down in a more expanded manner makes it more explicit:
\begin{align} \Bbb P(Y \in [a,b]) &= \Bbb P\Big(X \in g^{-1}([a,b])\Big) \\&= \Bbb P \Big(\cup_k \big(X \in g^{-1}([a,b]) \cap U_k\big)\Big) \\&= \sum_k \Bbb P(X \in g^{-1}([a,b]) \cap U_k) \\&= \sum_k \Bbb P\Big(X \in g|_k^{-1}([a,b]) \Big) \\&= \sum_k \Bbb P\Big(\underbrace{g|_k(X)}_{Y_k} \in [a,b] \Big) \\&= \sum_k \int_{[a,b]}f_{Y_k}(y)\,dy \\&= \sum_k \int_{[a,b]} f_X (g^{-1}|_k(y)) \left|\frac{\mathrm{d}}{\mathrm{d}y} g^{-1}|_k(y)\right| \,dy \tag{1} \\&= \int_{[a,b]} \sum_k f_X (g^{-1}|_k(y)) \left|\frac{\mathrm{d}}{\mathrm{d}y} g^{-1}|_k(y)\right| \,dy \end{align}
Notice that at $(1)$ we used the previous result.