Showing the relationship between $B$ and $\Gamma$ without multivariate calculus I have already read Gamma and Beta function proof.
Here is my work. By definition, the Beta function is 
$$B(a, b) = \int_{0}^{1}x^{a-1}(1-x)^{b-1}\text{ d}x$$
Let $y = \dfrac{x}{1-x}$, the "odds ratio" of $x$. Then after some algebra, it can be shown that $x = \dfrac{y}{1+y}$, $\text{d}x = \dfrac{1}{(1+y)^2}\text{ d}y$, and $1-x=\dfrac{1}{1+y}$, yielding
$$B(a, b) = \int_{0}^{\infty}\dfrac{y^{a-1}}{(1+y)^{a+b}}\text{ d}y\text{.}$$
Now we have, for the Gamma function, 
$$\Gamma(a) = \int_{0}^{\infty}y^{a-1}e^{-y}\text{ d}y\text{.}$$
For $k > 0$, let $ku = y$, so that $k\text{ d}u = \text{d}y$, hence after some work, we obtain
$$\Gamma(a) = k^a\int_{0}^{\infty}u^{a-1}e^{-ku}\text{ d}u$$
from which we obtain
$$\dfrac{1}{k^a}=\dfrac{\displaystyle \int_{0}^{\infty}u^{a-1}e^{-ku}\text{ d}u}{\Gamma(a)}$$
Now, the linked question above says - if I'm interpreting it correctly - to notice
$$\dfrac{1}{(1+y)^{a+b}}=\dfrac{\displaystyle \int_{0}^{\infty}u^{a+b-1}e^{-(1+y)u}\text{ d}u}{\Gamma(a+b)}$$
and then substitute this into $B(a, b)$:
$$\begin{align}B(a, b) &= \dfrac{1}{\Gamma(a+b)}\int_{0}^{\infty}y^{a-1}\left[\int_{0}^{\infty}u^{a+b-1}e^{-(1+y)u}\text{ d}u \right]\text{ d}y \\
&= \dfrac{1}{\Gamma(a+b)}\int_{0}^{\infty}y^{a+b-1}\left[\dfrac{\Gamma(a+b)}{(1+y)^{a+b}} \right]\text{ d}y
\end{align}$$
... and I'm not sure how this helps me.
 A: You have a double integral at some point, so your approach is using multivariate calculus.
But I can outline an approach which really avoids it. For any $a,b>0$, it is pretty clear that $B(a,b)=B(b,a)$ and that both $B(\cdot,b), B(a,\cdot)$ are continuous, non-negative functions on $\mathbb{R}^+$. We have $B(a,1)=\frac{1}{a}$ and by integration by parts we also have
$$ \frac{B(a+1,b)}{B(a,b)}= \frac{a}{a+b}$$
so for a fixed value of $b\in\mathbb{R}^+$ the functions
$$ f(a) = B(a,b)\Gamma(a+b),\qquad g(a)=\Gamma(a)\Gamma(b)$$
fulfill $f(1)=g(1)=\Gamma(b)$ and $\frac{f(a+1)}{f(a)}=\frac{g(a+1)}{g(a)}=a$ for any $a\in\mathbb{R}^+$.
By the Bohr-Mollerup theorem $\Gamma$ is log-convex on $\mathbb{R}^+$.
By the Cauchy-Schwarz inequality $B(a,b)$ is midpoint-log-convex with respect to the $a$ variable:
$$ \int_{0}^{1}x^{a_1}(1-x)^b\frac{dx}{x(1-x)}\int_{0}^{1}x^{a_2}(1-x)^b\frac{dx}{x(1-x)} \geq \left(\int_{0}^{1}x^{\frac{a_1+a_2}{2}}(1-x)^b\frac{dx}{x(1-x)}\right)^2$$
and since $B(a,b)$ is non-negative and continuous, it is log-convex.
By the Bohr-Mollerup theorem again, $f(a)$ and $g(a)$ are the same function on $\mathbb{R}^+$, hence
$$ B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} $$
as wanted.
