The Euler-Bernoulli differential equation without load is:
$$ EI{\frac {\partial ^{4}w}{\partial x^{4}}}+\mu {\frac {\partial ^{2}w}{\partial t^{2}}}=0 \tag{1}$$
Using $w(x,t)=\phi(x)q(t)$ the solutions for a cantilevered beam are:
$$\phi_n(x) = A_n\left[(\cosh \beta _nx-\cos \beta _nx) + \frac {\cos \beta _{n}L+\cosh \beta _{n}L}{\sin \beta _{n}L+\sinh \beta _{n}L}(\sin \beta _nx-\sinh \beta _nx)\right]$$
$$q(t)=a\sin(\omega_n t)+b\cos(\omega_n t)$$
where $L\in \mathbb{R}^+$ and $\beta_n$ is one of the positive solutions ot the equation:
$$ 1+\cos(\beta L)\cosh(\beta L)=0 \tag{2}$$
Then $\phi_n(x)$ are the eigenfunctions of the differential operator $\mathcal{D}= EI{\frac {\partial ^{4}}{\partial x^{4}}}-\omega_n^2\mu$:
$$ \mathcal{D}\phi_n(x)=\left(EI\beta_n^4-\omega_n^2 \mu\right)\phi_n(x)=\lambda_n\phi_n(x)$$
How can I prove that these functions are orthogonal?
$$ \int_0^L\phi_n(x)\phi_m(x) dx = L A_n^2 \delta_{nm} $$
Doing that explicitly is not possible because all the $\beta_n$ are irrational and I'm limited by the numerical precision. I did some numerical integrations that suggested me the above result, but do you know if there is quick way to prove that for every $n$?
I can also see by direct calculations that $\phi_n(L)=2A_n(-1)^{n+1}$ for $n=1,2,3,4$ but how can I prove that for every $n$?
This is done by partial integration, four times. We have $$\begin{align} \phi_n'''' &= \beta_n^4\phi_n & \phi_n(0) &= 0 & \phi_n'''(L) &= 0 \\&&\phi_n'(0) &= 0 & \phi_n''(L) &= 0 \end{align}$$ Therefore (I write $\phi_n$ instead of $\phi_n(x)$ for brevity): $$\begin{align} \beta_n^4\int_0^L \phi_n\,\phi_m\,\mathrm{d}x &= \int_0^L \phi_n''''\,\phi_m\,\mathrm{d}x \\&= \underbrace{\left[\phi_n'''\,\phi_m\right]_0^L}_0 - \underbrace{\left[\phi_n''\,\phi_m'\right]_0^L}_0 + \underbrace{\left[\phi_n'\,\phi_m''\right]_0^L}_0 - \underbrace{\left[\phi_n\,\phi_m'''\right]_0^L}_0 + \int_0^L \phi_n\,\phi_m''''\,\mathrm{d}x \\&= \beta_m^4\int_0^L \phi_n\,\phi_m\,\mathrm{d}x \end{align}$$ So for $\beta_n\neq\beta_m$, we must have $\int_0^L \phi_n\,\phi_m\,\mathrm{d}x = 0$.
The case $m=n$ is mere calculation of trigonometric integrals. I would (and should!) leave that to you. However, there are some ways to speed things up, and these deserve to be mentioned. So I make an exception here. A rare one.
For the $m=n$ case, it is best to become more familiar with the basis functions used for solving initial value problems with ODE $y''''=y$. Let $$\begin{align} R_0(z) &= \frac{1}{2}(\cosh z + \cos z) = R_1'(z) & R_2(z) &= \frac{1}{2}(\cosh z - \cos z) = R_3'(z) \\ R_1(z) &= \frac{1}{2}(\sinh z + \sin z) = R_2'(z) & R_3(z) &= \frac{1}{2}(\sinh z - \sin z) = R_0'(z) \end{align}$$ Then the trigonometric identities $\cos^2+\sin^2=1$ and $\cosh^2-\sinh^2=1$ turn into $$\begin{align} R_0^2 - 2\,R_1 R_3 + R_2^2 &= 1 \\ R_1^2 - 2\,R_0 R_2 + R_3^2 &= 0 \end{align}$$ To prove these, it is easiest to verify that the equations are fulfilled for $z=0$ and that differentiating the left hand sides yields zero.
Some integrals. The proof goes as before: Check at $z=0$ and differentiate. $$\begin{align} \int_0^z R_2^2(\zeta)\,\mathrm{d}\zeta &= \frac{1}{4}\Bigl(3\,R_2(z)\,R_3(z) - R_0(z)\,R_1(z) + z\Bigr) \\ \int_0^z R_2(\zeta)\,R_3(\zeta)\,\mathrm{d}\zeta &= \frac{1}{2}\,R_3^2(z) \\ \int_0^z R_3^2(\zeta)\,\mathrm{d}\zeta &= \frac{1}{4}\Bigl(3\,R_3(z)\,R_0(z) - R_1(z)\,R_2(z)\Bigr) \end{align}$$ These could be simplified further, but that would take doubled arguments like $(2z)$, and we would have to transform those back later. So let us stick to the above expressions.
For brevity, I will omit the subscript $n$ in $\phi_n,A_n,\beta_n$. We have $$\phi(x) = 2A\Bigl(R_2(\beta x) - B\,R_3(\beta x)\Bigr) \qquad B = \frac{R_0(\beta L)}{R_1(\beta L)} = \frac{R_3(\beta L)}{R_0(\beta L)}$$ The equations for $B$ follow from the boundary conditions $\phi''(L) = 0$ and $\phi'''(L) = 0$. These imply $$\underbrace{R_0^2(z) - R_1(z)\,R_3(z)} _{\frac{1}{2}\bigl(1 + \cosh(z)\cos(z)\bigr)} = 0 \quad\text{at}\quad z=\beta L$$ and thus lead to the characteristic equation $(2)$ mentioned in the question.
Integrating now yields $$\begin{align} \frac{1}{A^2}\int_0^L \phi^2(x)\,\mathrm{d}x &= 4\int_0^L \Bigl(R_2(\beta x) - B\,R_3(\beta x)\Bigr)^2\,\mathrm{d}x \\ &= 4\int_0^L R_2^2(\beta x)\,\mathrm{d}x - 8 B \int_0^L R_2(\beta x)\,R_3(\beta x)\,\mathrm{d}x + 4 B^2\int_0^L R_3^2(\beta x)\,\mathrm{d}x \\ &= \frac{1}{\beta}\Bigl(3\,R_2(\beta L)\,R_3(\beta L) - R_0(\beta L)\,R_1(\beta L) + \beta L\Bigr) - \frac{4 B}{\beta} R_3^2(\beta L) \\&\quad{} + \frac{B^2}{\beta}\Bigl(3\,R_3(\beta L)\,R_0(\beta L) - R_1(\beta L)\,R_2(\beta L)\Bigr) \end{align}$$ Move $B$ or $B^{-1}$ into the big parentheses, always using the form that cancels denominators: $$\begin{align} \frac{1}{A^2}\int_0^L \phi^2(x)\,\mathrm{d}x &= L + \frac{B}{\beta}\Bigl(3\,R_0(\beta L)\,R_2(\beta L) - R_1^2(\beta L)\Bigr) - \frac{4B}{\beta} R_3^2(\beta L) \\&\qquad{} + \frac{B}{\beta}\Bigl(3\,R_3^2(\beta L) - R_0(\beta L)\,R_2(\beta L)\Bigr) \\ &= L + \frac{B}{\beta}\underbrace{\Bigl( 2\,R_0(\beta L)\,R_2(\beta L) - R_1^2(\beta L) - R_3^2(\beta L)\Bigr)} _{0\quad\text{(trigonometric!)}} \end{align}$$ which proves that $\int_0^L \phi^2(x)\,\mathrm{d}x = A^2 L$.
The moral of the story: When expressions get complicated:
Here is a useful exercise: Apply other zero-energy boundary conditions, verify that the orthogonality remains intact, and calculate the $m=n$ integral for those boundary conditions.