How can one prove the equation for a system of lenses? Let's say we have a system of 2 lenses which are placed in a non-negligible distance $L$ from eachother. 
How can I prove without using geometrical construction and principal planes that
$$\frac 1f = \frac {1}{f_1} + \frac {1}{f_2} -\frac{L}{f_1 f_2} $$
describes the system's focal length?

For example, with $ L \lt \lt f_1f_2 $ :
$\dfrac{1}{f_1} = \dfrac{1}{o_1} + \dfrac{1}{i_1}$ where $o_1$ is the object distance and $i_1$ the image distance of the first lens. 
And of course, $\dfrac{1}{f_2} = \dfrac{1}{o_2} + \dfrac{1}{i_2}$ .
Now you can say the image of the first lense is the object distance of the second length (with paying attention to conventions what is negative and what positive):
$o_2 = -i_1$ 
Using this relation one receives $\dfrac{1}{o_1} + \dfrac{1}{i_2} = \dfrac{1}{f_1} + \dfrac{1}{f_2}  = \dfrac{1}{f}$ 

I already tried saying $o_2 = -i_1 + L$ but this didn't lead anywhere.
 A: I take it you mean you wish to look at this problem without ray constructions. To do it without "geometrical" constructions is impossible because this is fundamentally a problem about geometry! That is, we want to know where and how fields converge, where they are plane and so forth, so I'd argue you're using geometrical arguments even if you analyze this problem by solving Maxwell's equations.
One way to do this is to think of the fields as spherical waves (with a plane wave as a zero curvature case thereof). In this picture, the lenses are operators that act on and transform the curvatures of the input wavefront and output a different curvature spherical wave. Likewise, propagation through freespace is also such a transformer: a converging spherical wavefront of radius $R$ becomes one of radius $R-d$ if it propagates a distance $d$ towards its focus.
To make this method work, you need to think of the lens as a phase mask that adds / subtracts curvature to / from spherical wavefronts. The details of what this statement means and why it is true (from a scalar wave theory standpoint) are given in my answer here. 
I shall now think of the field propagating left to right and use positive curvatures for wavefronts converging to a focus and negative for those diverging from one.
So we begin with a spherical wave diverging from a point $d_i$ before the first lens. We therefore have a spherical wave of curvature $-1/d_i$ input to the first lens. The latter adds curvature $1/f_1$, so the its output is a converging wave, curvature $f_1^{-1}-d_i^{-1}$ (this could still be negative, thus diverging, but properly signed quantities take care of themselves. So the output radius of curvature is $(f_1^{-1}-d_i^{-1})^{-1}$. The spherical wavefront now propagates a distance $L$, so its radius is $(f_1^{-1}-d_i^{-1})^{-1}-L$ when it reaches the second lens, i.e. its curvature is $\left((f_1^{-1}-d_i^{-1})^{-1}-L\right)^{-1}$. The second lens adds curvature $1/f_2$. So the curvature of the output from the second lens is simply:
$$\kappa_2 = \frac{1}{f_d}=\left((f_1^{-1}-d_i^{-1})^{-1}-L\right)^{-1} + f_2^{-1}$$
Here $f_d$ is the distance beyond the second lens that the field must travel to reach its focus. The focal length is derived from this distance $f_d$ as $d_i\to \infty$ (which we could also have gotten from assuming a plane wave input), namely the quantity $f_d=((f_1-L)^{-1} + f_2^{-1})^{-1}$.
You can do analogously for fields propagating right to left and find that a plane wave input to the right hand side converges to the left of the lens with focal length $f_1$ at a distance $((f_2-L)^{-1} + f_1^{-1})^{-1}$ from that lens.
Now, even though you now know where the fields converge, these quantities are not the same as the focal length: to derive the formula you have cited, you need to know where the output principal plane is: the focal length is defined from this plane. So unfortunately here is one point where you need to use the formula to locate the principal planes, you can't avoid that. You may be helped by my explanation of the principal planes here.
A: Use ray transfer matrices. The system matrix is given by
$$
S=\begin{pmatrix}
1 & 0 \\
-\frac{1}{f_2} & 1 \\
\end{pmatrix}
\begin{pmatrix}
        1 & L \\
        0 & 1 \\
        \end{pmatrix}
\begin{pmatrix}
        1 & 0 \\
        -\frac{1}{f_1} & 1 \\
        \end{pmatrix}
$$
The back focal length $BFL$ is determined by the point where all the rays in an incoming collimated beam meet after the second lens:
$$
\begin{pmatrix}
        1 & BFL \\
        0 & 1 \\
        \end{pmatrix}
S
\begin{pmatrix}
        z \\
        0  \\
        \end{pmatrix}
=
\begin{pmatrix}
        0 \\
        \text{X}  \\
        \end{pmatrix}
$$
for all $z$. Solve for $BFL$ and then find $f$ by adding the distance of the second lens to the origin of coordinates. The origin is not arbitrary, it is chosen to a specific point on purpose in order to give a neat symmetric expression.
