Here is the proof$^\ast$ (and answer to my question) which I promised.
The functions $f_1$ and $f_2$ form the set of affine functions and we can define
\begin{align}
f_1(x) = \frac{1}{\sqrt{2}}, \quad f_2(x) = \sqrt{\frac{3}{2}} x.
\end{align}
All subsequent functions are given (in correct order and normalized) by
\begin{align}
\sqrt{1+\frac{\cos(\nu_{1+2i})^2}{\cosh(\nu_{1+2i})^2}} f_{1+2i} & = \cos(\nu_{1+2i}x)+\frac{\cos(\nu_{1+2i})}{\cosh(\nu_{1+2i})} \cosh(\nu_{1+2i}x) \\
\sqrt{1-\frac{\sin(\nu_{2+2i})^2}{\sinh(\nu_{2+2i})^2}} f_{2+2i} & = \sin(\nu_{2+2i}x)+\frac{\sin(\nu_{2+2i})}{\sinh(\nu_{2+2i})} \sinh(\nu_{2+2i}x), \quad i \in \mathbb{N}
\end{align}
where $\nu_{1+2i}$ is the $i$-th (positive) root of $\cos(\nu)\sinh(\nu)+\cosh(\nu)\sin(\nu)$ and $\nu_{2+2i}$ is the $i$-th (positive) root of $\cos(\nu)\sinh(\nu)-\cosh(\nu)\sin(\nu)$. Further more, we have
\begin{align}
\int_{-1}^1 f_i'' f_j'' = \delta_{ij} \nu_i^4.
\end{align}
$^\ast$: At one point we assume $f_k$ to be sufficiently smooth. I guess this follows easily within a standard argument using calculus of variations, but every derivation I found so far, just makes the same assumption.
Proof:
Using Lagrange Multipliers (this part is thanks to Dirk) in order to obtain $f_k$, $k > 3$, we form the Lagrangian
\begin{align}
\mathcal{L}^{(k)}(g,\lambda^{(k)},\mu^{(k)}) & = \int_{-1}^1 |g''|^2 + \sum_{i = 1}^{k-1} \lambda_i^{(k)} \int_{-1}^1 g f_i
+ \mu^{(k)} \left( \int_{-1}^1 |g|^2 - 1 \right).
\end{align}
The derivatives with respect to $\lambda_1,\ldots,\lambda_{k-1}$ give the orthogonality conditions,
whereas $\mu$ gives the normalization. Variation of the function $g$ provides
\begin{align}
2 \int_{-1}^1 g'' h'' + \sum_{i = 1}^{k-1} \lambda_i^{(k)} \int_{-1}^1 f_i h
+ 2 \mu^{(k)} \int_{-1}^1 g h = 0.
\end{align}
Using partial integration twice (at this points we actually assume $f$ to be sufficiently smooth)
we get
\begin{align}
& 2 \int_{-1}^1 g^{(4)} h + 2 g^{(2)} h|_{-1}^1 - 2 g^{(3)} h|_{-1}^1 \\
& + \sum_{i = 1}^{k-1} \lambda_i^{(k)} \int_{-1}^1 f_i h
+ 2 \mu^{(k)} \int_{-1}^1 g h = 0.
\end{align}
Since $h \in C^\infty$ is arbitrary, we get the pointwise identity
\begin{align}
\mathcal{P}^{(k)}(g) := 2g^{(4)} + \sum_{i = 1}^{k-1} \lambda_i^{(k)}f_i
+ 2 \mu^{(k)} g = 0
\end{align}
as well as the boundary conditions
\begin{align}
g^{(2)}(\pm 1) = 0, \quad g^{(3)}(\pm 1) = 0.
\end{align}
We inductively proof now that $\lambda^{(k)} \equiv 0$, $k \geq 3$. For $j < k$, by orthonormality conditions,
we can simplify
\begin{align}
0 = \int_{-1}^1 \mathcal{P}^{(k)}(g) f_j & =
2 \int_{-1}^1 g^{(4)} f_j + \sum_{i = 1}^{k-1} \lambda_i^{(k)} \int_{-1}^1 f_i f_j
+ 2 \int_{-1}^1 \mu^{(k)} g f_j \\
& =
2 \int_{-1}^1 g^{(4)} f_j + \lambda_j^{(k)}
\end{align}
Using the boundary conditions for both $f$ and $g$, this is equivalent to
\begin{align}
0 = \int_{-1}^1 g f_j^{(4)} + \lambda_j^{(k)}
\end{align}
In the case of $k = 3$, both $f_1^{(4)} \equiv f_2^{(4)} \equiv 0$.
Hence, $\lambda_j^{(3)} = 0$, $j = 1,2$. We remain with
\begin{align}
g^{(4)} + \mu^{(k)} g = 0
\end{align}
In the $k$-th inductive step, we can hence assume that $f_j^{(4)} = -\mu^{(j)} f_j$.
Again, by orthogonality, we conclude $\lambda_j^{(k)} = 0$ for (all) $k$. Thereby,
$f_k$ is the one normalized solution of the simple differential equation
\begin{align}
f_k^{(4)} + \mu^{(k)} f_k = 0, \quad f_k^{(2)}(\pm 1) = 0, \ f_k^{(3)}(\pm 1) = 0.
\end{align}
Searching for its solutions, $f_k$ can be assumed to be either symmetric or skew-symmetric.
If neither was the case, then due to the symmetry of the problem, $g(x) := f_k(-x)$ fulfills the same differential equation
(for the same $\mu^{(k)}$). Hence, $g + f_k$ and $g - f_k$ are
symmetric and skew-symmetric, respectively, fulfill the differential equation and span the same space
as $g$ and $f_k$. For simplicity, we substitute $\nu = \sqrt[4]{-\mu^{(k)}}$.
The unscaled symmetric and skew-symmetric solutions are given by
\begin{align}
\cos(\nu x) + a \cosh(\nu x), \quad \sin(\nu x) + b \sinh(\nu x),
\end{align}
respectively.
In order to satisfy the boundary conditions,
one calculates that
\begin{align}
\begin{pmatrix}
-\cos(\nu) & \cosh(\nu) \\ \sin(\nu) & \sinh(\nu)
\end{pmatrix}
\begin{pmatrix}
1 \\ a
\end{pmatrix}
= 0, \quad
\begin{pmatrix}
-\sin(\nu) & \sinh(\nu) \\ -\cos(\nu) & \cosh(\nu)
\end{pmatrix}
\begin{pmatrix}
1 \\ b
\end{pmatrix}
= 0.
\end{align}
The points at which either matrix becomes singular, are the roots $r_c$ and $r_s$ of
\begin{align}
t_c := \cos(\nu)\sinh(\nu)+\cosh(\nu)\sin(\nu), \quad t_s := \cos(\nu)\sinh(\nu)-\cosh(\nu)\sin(\nu),
\end{align}
respectively. Both have infinitely many (single) roots which can be approximated by $-\pi/4+i\pi$ and $\pi/4+i\pi$, respectively, since
for $k > 7$, these equal the true $i$-th roots $r_c^{(i)}$ and $r_s^{(i)}$ up to $20$ decimal places.
Earlier values are easily found numerically.
We therefor define $\nu_{1,2} = 0, \ \nu_{1+2i} = r^{(i)}_c$ and $\nu_{2+2i} = r^{(i)}_s$.
After some calculus, we obtain
\begin{align}
\sqrt{1+\frac{\cos(\nu_{1+2i})^2}{\cosh(\nu_{1+2i})^2}} f_{1+2i} & = \cos(\nu_{1+2i}x)+\frac{\cos(\nu_{1+2i})}{\cosh(\nu_{1+2i})} \cosh(\nu_{1+2i}x) \\
\sqrt{1-\frac{\sin(\nu_{2+2i})^2}{\sinh(\nu_{2+2i})^2}} f_{2+2i} & = \sin(\nu_{2+2i}x)+\frac{\sin(\nu_{2+2i})}{\sinh(\nu_{2+2i})} \sinh(\nu_{2+2i}x), \quad i \in \mathbb{N}
\end{align}
All these roots are distinct. The orthogonality conditions are easy to verify, considering
\begin{align}
\int_{-1}^1 f_i f_j = \frac{1}{{\nu_i}^4} \int_{-1}^1 f_i^{(4)} f_j = -\frac{1}{{\nu_i}^4} \int_{-1}^1 f_i f_j^{(4)}
= \frac{\nu_j^4}{{\nu_i}^4} \int_{-1}^1 f_i f_j
\end{align}
since either $\nu_j \neq \nu_i$ or $\nu_j = 0$, $j = 1,2$, for $j < i$.
One can further calculate
\begin{align}
\int_{-1}^1 f_i'' f_j'' = \delta_{ij} \nu_i^4.
\end{align}
The $f_k$ are hence in correct order and their second derivates are orthogonal as well.
