Let $f$ be a continuous real-valued bounded function defined on $[0,+\infty[$.
We want to compute $\mathbb{E}\big[ f(Z) \big]$. Given that $X$ and $Y$ are independent (and both follow exponential distributions), we have:
$$ \mathbb{E}\big[ f(Z) \big] = \int_{0}^{+\infty} \int_{0}^{+\infty} f\Big(\frac{x}{y+1}\Big) \lambda_1 e^{-\lambda_{1}x} \lambda_{2}e^{-\lambda_{2}y} \; dxdy. $$
Consider the change of variables $\Phi \, : \, [0,+\infty[^2 \, \rightarrow \, [0,+\infty[^2$ defined by:
$$ \forall (x,y) \in [0,+\infty[^2, \; \Phi\big( (x,y) \big) = \Big(\frac{x}{y+1}, y\Big) = (u,v) \in [0,+\infty[^2. $$
$\Phi$ is a $\mathcal{C}^{1}$ diffeomorphism from $[0,+\infty[^2$ to itself and:
$$ \Phi^{-1}\big( (u,v) \big) = \big( u(v+1), v \big) $$
$$ \mathrm{Jac}\big( \Phi^{-1}, (u,v) \big) = \begin{bmatrix} v+1 & u \\ 0 & 1 \end{bmatrix}. $$
It follows from the change of variable theorem that:
$$
\begin{align*}
\mathbb{E}\big[ f(Z) \big] & = {} \int_{0}^{+\infty} \int_{0}^{+\infty} f(u) \lambda_1\lambda_2 e^{-\lambda_{1}u(v+1)}e^{-\lambda_{2}v} (v+1) \, du dv. \\[2mm]
& = \int_{0}^{+\infty} f(u) \Bigg( \lambda_1 \lambda_2 e^{-\lambda_1 u} \int_{0}^{+\infty} (v+1) e^{-(\lambda_{1}u + \lambda_{2})v} \, dv \Bigg).
\end{align*}$$
Also, note that:
$$
\begin{align*}
\int_{0}^{+\infty} (v+1)e^{-(\lambda_{1}u + \lambda_{2})v} dv & = {} \Bigg[-\frac{e^{-av}(av + a + 1)}{a^2} \Bigg]_{v=0}^{v=+\infty} \\
& = \frac{a+1}{a^2}.
\end{align*}$$
with $a = \lambda_{1}u + \lambda_{2}$.
Therefore, the probability density function of $Z$ (with respect to the Lebesgue measure on $[0,+\infty[$) is given by:
$$ u \in [0,+\infty[ \, \mapsto \; \lambda_{1}\lambda_{2}e^{-\lambda_{1}u} \frac{\lambda_{1}u + \lambda_{2} + 1}{(\lambda_{1}u + \lambda_{2})^2}. $$
Assuming you mean $\lambda$ is the rate parameter here (i.e. Exponential with mean $1/\lambda$).
First of all, recheck your density of $S$. The correct density as mentioned in comments is $$f_S(s)=\lambda^2se^{-\lambda s}\mathbf1_{s>0}$$
It is easy to verify that the density of $T$ is $$f_T(t)=\mathbf1_{0<t<1}$$
You can find the joint distribution function of $(S,T)$ as follows:
For $s>0$ and $0<t<1$,
\begin{align}
P(S\le s,T\le t)&=P\left(X+Y\le s,\frac{X}{X+Y}\le t\right)
\\&=\iint_D f_{X,Y}(x,y)\,dx\,dy\quad\quad,\text{ where }D=\{(x,y):x+y\le s,x/(x+y)\le t\}
\\&=\lambda^2\iint_D e^{-\lambda(x+y)}\mathbf1_{x,y>0}\,dx\,dy
\end{align}
Change variables $(x,y)\to(u,v)$ such that $$u=x+y\quad,\quad v=\frac{x}{x+y}$$
This implies $$x=uv\quad,\quad y=u(1-v)$$
Clearly, $$x,y>0\implies u>0\,,\,0<v<1$$
And $$dx\,dy=u\,du\,dv$$
So again for $s>0\,,\,0<t<1$,
\begin{align}
P(S\le s,T\le t)&=\lambda^2\iint_R ue^{-\lambda u}\mathbf1_{u>0,0<v<1}\,du\,dv\qquad,\text{ where }R=\{(u,v):u\le s,v\le t\}
\\&=\left(\int_0^s \lambda^2 ue^{-\lambda u}\,du\right) \left(\int_0^t \,dv\right)
\\\\&=P(S\le s)\,P(T\le t)
\end{align}
This proves the independence of $S$ and $T$, with $S$ a Gamma variable and $T$ a $U(0,1)$ variable.
And from the joint distribution function, it is readily seen (without differentiating) that the joint density of $(S,T)$ is $$f_{S,T}(s,t)=\lambda^2 se^{-\lambda s}\mathbf1_{s>0,0<t<1}$$
Note that the change of variables isn't really necessary if you are comfortable with the first form of the double integral.
Best Answer
Let $W:=(Y,Z)$, then as $Y$ and $Z$ are independent $f_W=f_Y\cdot f_Z$, therefore for $A_c:=\{(s,t)\in \mathbb{R}^2 :s/t\leqslant c\}$ we have that
$$ \begin{align*} \Pr [Z/Y\leqslant c]&=\Pr [W\in A_c]\\&=\int_{\mathbb{R}^2}\mathbf{1}_{A_c}(s,t)f_Y(s)f_Z(t)\mathop{}\!d(s,t)\\ &=\frac1{2\pi }\left(\int_{(0,\infty )}e^{-t^2/2}\int_{(-\infty ,tc]}e^{-s^2/2}\mathop{}\!d s \mathop{}\!d t+\int_{(-\infty ,0)}e^{-t^2/2}\int_{[tc,\infty )}e^{-s^2/2}\mathop{}\!d s \mathop{}\!d t\right) \end{align*} $$ Differentiating respect to $c$ gives $$ f_{Z/Y}(c)=\frac1{2\pi}\left(\int_{(0,\infty )}te^{-\frac1{2}(1+c)^2t^2}\mathop{}\!d t-\int_{(-\infty ,0)}te^{-\frac1{2}(1+c)^2t^2}\mathop{}\!d t\right)\\ =\frac1{\pi }\int_{(0,\infty )}te^{-\frac1{2}(1+c)^2t^2}\mathop{}\!d t=\frac1{\pi(1+c^2)} $$ ∎
EDIT: alternatively, using the change to polar coordinates $(t,s)=(r\cos \alpha ,r\sin \alpha )$ we have that the condition $s/t\leqslant c$ is equivalent to $\tan \alpha \leqslant c \Rightarrow -\pi/2<\alpha \leqslant \arctan c$, and as the tangent have two full periods on $(-\pi,\pi)$ then
$$ \int_{\mathbb{R}^2}\mathbf{1}_{A_c}(s,t)f_Y(s)f_Z(t)\mathop{}\!d(s,t)=\frac1{\pi }\int_{[0,\infty )\times (-\pi/2,\arctan c]}re^{-r^2/2}\mathop{}\!d (r,\alpha ) =\frac1{\pi}(\arctan c+\frac{\pi}{2}) $$
where the result follows again by differentiation.