(*** w(0) < d0  with P > 2*Sqrt[Dk]  **)
Print["y(x)"]
y=b1*Cos[alp1*x]+b2*Cos[bet1*x]
Print["Differential Equation"]
diffeq=FullSimplify[dd*D[y,{x,4}]+ph*D[y,{x,2}]+k*y]
alp1=Sqrt[(ph+Sqrt[ph^2-4*dd*k])/(2*dd)];
bet1=Sqrt[(ph-Sqrt[ph^2-4*dd*k])/(2*dd)];
FullSimplify[diffeq, Assumptions -> {k>0,dd>0}]
alp1=.
bet1=.
Print["delta*L"]
deltaL=Simplify[ell*Integrate[(D[y,x])^2,{x,0,ell}]/2 /. {alp1->alp1L/ell,bet1->bet1L/ell,c->cL*ell}] 
Print["Two BCs"]
bc2=y /. x->ell;
bc2=bc2 /. {alp1->alp1L/ell,bet1->bet1L/ell};
bc2=Simplify[bc2 /. c-> cL*ell]
bc3=D[y,x] /. x->ell;
bc3=ell*bc3 /. {alp1->alp1L/ell,bet1->bet1L/ell};
bc3=Simplify[bc3 /. c-> cL*ell]
a11=bc2 /. {b1->1,b2->0};
a12=bc2 /. {b1->0,b2->1};
a21=bc3 /. {b1->1,b2->0};
a22=bc3 /. {b1->0,b2->1};
Print["Determinant"]
det=a11*a22-a21*a12
Print["b1/b2"]
b12=-a12/a11
b1=b12*b2;
Print["det=0 Condition"]
bc1=b1+b2-w0;
sol2=Solve[bc1==0,b2];
b2=b2 /. sol2;
b2=Part[b2,1];
eq1=det/. {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]]}
Print["Moment BC at End"]
bc4a=ell^2*D[y,{x,2}] /. x->ell;
bc4a=bc4a /. {alp1->alp1L/ell,bet1->bet1L/ell};
bc4a=Simplify[bc4a/d0 /. {c-> cL*ell}];
bc4b=Simplify[6*deltaL/(d0^2*hd)];
bc4=Simplify[bc4a-bc4b];
bc4=Simplify[bc4 /. w0->w0d0*d0];
bc4=Simplify[bc4/w0d0]
eq2=bc4/. {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]]}
Save["tempfile1.mth",{eq1,eq2,deltaL}]

<<"tempfile1.mth"; 
m=9.699
hdd=25
phLDG=12.0;
kLDG=2.6^4;
wstart=0.001;
iiimax=11;
winc=(1-wstart)/(iiimax-1);
w0d0s=wstart-winc;
Print[" w0/d0   phL^2/2D  ph/2(kD)^(1/2)  L(k/D)^(1/4)  delph/2(kD)^(1/2)  Energy Frac"]
For[iii=0,iii<iiimax,iii++,w0d0s=w0d0s+winc;
eq1c=eq1 /. {w0d0->w0d0s};
eq2c=eq2 /. {w0d0->w0d0s,hd->hdd};
sol2=FindRoot[{eq1c==0,eq2c==0},{{phLD,phLDG},{kLD,kLDG}}];
pLs=phLD /. sol2; phLDG=pLs;
kLs=kLD /. sol2; kLDG=kLs;
delphkD=6*deltaL/(hd^2*Sqrt[kLs]) /. {w0 -> w0d0s,hd->hdd,alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],phLD->pLs,kLD->kLs};
delphkD=delphkD /. {phLD->pLs,kLD->kLs};
enfrac=(1/3)*delphkD^2*hd^2/m /. hd->hdd;
Print[w0d0s,",  ",pLs,",   ",pLs/Sqrt[kLs]",     ",kLs^0.25,",     ",delphkD,",     ",enfrac,";"]]