Solution to invariance equation when deriving Lorentz transformation Requiring that the spacetime interval (1 spatial dimension) between the origin and an event $(t,\, x)$ stays constant under a transformation between reference frames: 
$$c^2t^2-x^2= c^2t'^2- x'^2$$
we want to find a solution to that equation, i.e. express $(t,\, x)$ in terms of $(t',\,x')$. In Landau/Lifshitz II it is stated that the general solution is
$$x= x' \mathrm{cosh} (\psi)+ ct' \mathrm{sinh} (\psi),\quad ct= x' \mathrm{sinh}(\psi) + ct' \mathrm{cosh}(\psi)$$
with $\psi$ being the rotating angle in the $tx$ - plane.
It is clear to me that the equation is then satisfied, but how would one find this solution starting from the given equation and why is it the most general one?

The different types of Lorentz transformations can be categorized using group theory, but in this case, I think one can convince themselves with a bit of algebra.
Let's write our transformation in matrix notation. Representing our event as a column vector, $\begin{pmatrix} ct \\ x \end{pmatrix}$, the transformation you wrote above is given as
$\begin{pmatrix}
\cosh (\psi) & \sinh(\psi) \\
\sinh (\psi) & \cosh(\psi)
\end{pmatrix}$
since 
$\begin{pmatrix}
\cosh (\psi) & \sinh(\psi) \\
\sinh (\psi) & \cosh(\psi)
\end{pmatrix}
\begin{pmatrix} ct \\ x \end{pmatrix} = 
\begin{pmatrix}
ct \cosh(\psi) + x \sinh(\psi)  \\ ct \sinh(\psi) + x \cosh(\psi)
\end{pmatrix}\,.$
More generally, transformations can be represented by 2 x 2 matrices,
$\begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} ct \\ x \end{pmatrix}\,.$
The spacetime interval is given by
$\begin{pmatrix} ct & x \end{pmatrix}
\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}
\begin{pmatrix} ct \\ x \end{pmatrix} = c^2 t^2 - x^2\,.$
Let's act on this value by a general transformation
$\begin{pmatrix} ct & x \end{pmatrix}
\begin{pmatrix} a & c \\ b & d\end{pmatrix}
\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}
\begin{pmatrix} a & b \\ c & d\end{pmatrix}
\begin{pmatrix} ct \\ x \end{pmatrix}\,.$
Note that the transformation acts twice - once for each copy of the event. Demanding that this value stay the same is identical to demanding 
$\begin{pmatrix} a & c \\ b & d\end{pmatrix}
\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}
\begin{pmatrix} a & b \\ c & d\end{pmatrix}
= \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\,.$
Doing out the matrix multiplication, we get the constraints $a^2 - c^2 = 1$, $b^2 - d^2 = -1$, and $ac-bd=0$. Assuming we are rotating in the $tx$ plane, $b$ and $c$ cannot equal zero, and we can rule out simple solutions like $a=\pm 1$, $d = \pm 1$. As far as I can tell, hyperbolic functions are the only functions that satisfy those above constraints.
Note that through this formulation, we've essentially reduced our physics question of Lorentz transformations to a mathematical question of matrices - specifically, what sort of matrices satisfy the last matrix equation. Matrices that do this are part of the indefinite orthogonal group $O(p,q)$ In this case, we have one space and one time component, so our group is $O(1,1)$ - I found a more mathematical discussion of how to get all the members of $O(1,1)$ here, in case the above wasn't sufficiently convincing. 
