I once did a drive to derive a symbolic solution (without trigonometry functions) for myself (mostly relying on Wolfram Alpha to do the heavy lifting) using the same constraints as @a-doering ($f(\pm1)=\pm1$ and $f'''(\pm1)=0$). I arrived at the fairly – all things considered – nice looking coefficients:
$\begin{align}
f(x) &= a \tanh(bx)
\\
a &= \sqrt{3} &&\approx 1.732050808 &&\approx 1.7159
\\
b &= \frac{-\ln(2 - \sqrt{3})}{2} &&\approx 0.658478948 &&\approx \frac{2}{3}
\end{align}$
Unfortunately, I do not remember the steps I took to get there.
Here's an interactive graph for anyone interested in playing around with it: https://www.desmos.com/calculator/tf4udjl8cn Original in red, mine in green. I doubt there would be any real difference in training outcome between them.
Addendum
I went back and rederived it, so here goes.
Start by deriving $b$ by setting $x=1$ and calculating $\frac{d^3}{dx^3} \tanh(b) = 0$ (the $a$-factor has no impact on the third derivative). WolframAlpha gives $b = \pm\frac{1}{2}\cosh^{-1}(2) = \pm\frac{1}{2}\ln(2+\sqrt{3})$.
Since setting $b=\frac{1}{2}\ln(2+\sqrt{3})$ guarantees that the third-derivative-condition will be met when $x=1$, we can lock both $x$ and $b$, and scale the output to fulfill the $f(1) = 1$ condition. Trusting in WolframAlpha, we have:
\begin{align}
a \tanh\left(\frac{\ln(2+\sqrt{3})}{2}\right) &= 1
\\
a &= \frac{1}{\tanh\left(\frac{\ln(2+\sqrt{3})}{2}\right)}
\\
&= \sqrt{3}
\end{align}
To see the equivalence with my original answer, note that $\ln(2+\sqrt{3}) = -\ln(2-\sqrt{3})$. The full function thus becomes:
\begin{align}
f(x) &= \sqrt{3} \tanh\left( \frac{\ln(2 + \sqrt{3})}{2} x \right)
\end{align}
