@Axel Vogt thanks for solution 

Since you say its numerically very ugly, maybe there is a proper way to deal with my problem. Basically I have to find omega by solving the following matrix determinant:

M := Matrix(4, 4, {(1, 1) = 0, (1, 2) = -1, (1, 3) = 0, (1, 4) = 1, (2, 1) = -EI*beta^3, (2, 2) = -m*omega^2, (2, 3) = EI*beta^3, (2, 4) = -m*omega^2, (3, 1) = EI*sin(147*beta)*beta+k_r*cos(147*beta)+I*c_r*omega*cos(147*beta), (3, 2) = EI*cos(147*beta)*beta-k_r*sin(147*beta)-I*c_r*omega*sin(147*beta), (3, 3) = -EI*sinh(147*beta)*beta+k_r*cosh(147*beta)+I*c_r*omega*cosh(147*beta), (3, 4) = -EI*cosh(147*beta)*beta+k_r*sinh(147*beta)+I*c_r*omega*sinh(147*beta), (4, 1) = -EI*cos(147*beta)*beta^3+sin(147*beta)*k_h, (4, 2) = EI*sin(147*beta)*beta^3+cos(147*beta)*k_h, (4, 3) = EI*cosh(147*beta)*beta^3+sinh(147*beta)*k_h, (4, 4) = EI*sinh(147*beta)*beta^3+cosh(147*beta)*k_h}):

The remaining values are:

beta=((5000/10^12)*(omega^2))^(1/4),k_r=3.33*10^10,k_h=1.62*10^9,c_r=3.14*10^9,m=350000,L=147,EI=10^12:

So maybe there is more decent ways to solve it numerically, or even a reasonable analytical expression?

@Preben Alsholm I dont get why D(PHI)(20) doesnt give me a value?

@Preben Alsholm mode shapes and eigen frequencies look right. Can you also express a function of phi(x) for omega1?

@Preben Alsholm k_r=k_h/2. As for eigenvalue, I need only first positive value,that will correspond first natural frequency of the beam.

@Preben Alsholm Indeed thats an eigenvalue problem. The difficulty of the problem is that two unknowns must be found- first mode shape (phi(x) and first eigenfrequency (omega/lambda). I assume some search routine has to be set up to find the initial unknown boundary condition and natural frequency so that it would match exact boundary conditions. For instance phi(L)'''=phi(L)'' then problem has to be worked out for several times modifying incorect values thus reducing the error. So some way of measuring the error has to be instituted. 

@Carl Love try this one. problem.mw. Let me know if it doesn't work

Wind.mw

@Carl Love this is the one that doesn't work: Wind.mw

Wind.mw                                      

Here is my complete worksheet. My orginal problem is to solve diff equations with 8 interface and boundary conditions. I wasn't able to do it with dsolve, so I did it manualy by substituting general soluctions to interface/boundary cond. It seems to be work fine as longs as my V1 and V2 are assumed to be 0. Most perfect scenario would be if I had solution with V1 and V2 being as a functions of z, or at least a constant values.