Integration problem : difference between two draws from an hyperbolic secant probability distribution I would like to find the distribution of the difference between two independent draws from a random variable which follows a secant hyperbolic distribution whose density function is :
$f(x) = \frac{1}{2}\operatorname{sech}(\frac{\pi}{2}x)$
where $\operatorname{sech}(x)=\cosh(x)^{-1}$
$\operatorname{sech}$ can be integrated (cf Wikipedia) :
$\int \operatorname{sech}(ax) = a^{-1} \arctan(\sinh(ax)) +C$
To get my distribution, I need to integrate :
$$g(x)=\int \operatorname{sech}(t-\frac{x}{2})\operatorname{sech}(t+\frac{x}{2})dt$$
$$g(x)=\int \frac{4}{e^{2t}+e^{-2t}+e^x+e^{-x}}dt$$
So, to put it simplier, I need to integrate :
$$h(x)=\frac{1}{e^x+e^{-x}+C}$$
I don't know how to go further.
We can notice that, when C=0, it is the same problem as integrating $\operatorname{sech}$ (Unfortunaltely, I have only the result, not the demonstration that could probably be helpful)
 A: Using Mathematica one can obtain $g(x)$:
f[x_] := 1/2 Sech[(x π)/2]
g[z_] := Evaluate[Integrate[f[x1] f[x1 + z], {x1, -∞, ∞}, Assumptions -> z ∈ Reals]]
g[z]
(* 1/2 z Csch[(π z)/2] *)

As a partial check take a random sample for $x_1$ and $x_2$, take the difference, and compare the resulting histogram to the probability density function:
dist = ProbabilityDistribution[(1/2) Sech[(x π)/2], {x, -∞, ∞}];
x1 = RandomVariate[dist, 1000000];
x2 = RandomVariate[dist, 1000000];
zz = x1 - x2;
Show[Histogram[zz, Automatic, "PDF"],
  Plot[g[z], {z, -5, 5}]]


** Steps to obtain the integral **
Using the Rubi package in Mathematica the first steps in solving the integral are as follows:
Steps[Int[f[x1] f[x1 + z], x1]]


This results in
-((Csch[(π z)/2] Log[Cosh[(π x1)/2]])/(2 π)) + (Csch[(π z)/2] Log[Cosh[(π x1)/2 + (π z)/2]])/(2 π)

for the anti-derivative.  The difference in the anti-derivative for $x_1=\infty$ and $x_1=-\infty$ is
-((Csch[(π z)/2] Log[E^(-((π z)/2))])/(2 π)) + (Csch[(π z)/2] Log[E^((π z)/2)])/(2 π)

This is simplified in Mathematica with
-((Csch[(π z)/2] Log[E^(-((π z)/2))])/(2 π)) + (Csch[(π z)/2] Log[E^((π z)/2)])/(2 π) /. Log[Exp[x_]] -> x

resulting in the answer:
1/2 z Csch[(π z)/2]

