(*** w(0)> d0  and w(0) < d1=m*d0  **)
rhs=-d0*k*m/(m-1);
kp=k/(m-1);
Print["y1(x)"]
y1=b7*Cos[alp1*x]+b8*Sin[alp1*x]+b9*Cos[bet1*x]+b10*Sin[bet1*x]
diffeq1=FullSimplify[dd*D[y1,{x,4}]+ph*D[y1,{x,2}]+k*y1];
alp1=Sqrt[(ph+Sqrt[ph^2-4*dd*k])/(2*dd)];
bet1=Sqrt[(ph-Sqrt[ph^2-4*dd*k])/(2*dd)];
Print["diffeq1"]
FullSimplify[diffeq1, Assumptions -> {k>0,dd>0}]
alp1=.
bet1=.
Print["y2(x)"]
y2=b3*Cos[alp2*x]+b4*Cosh[bet2*x]-rhs/kp
diffeq2=FullSimplify[dd*D[y2,{x,4}]+ph*D[y2,{x,2}]-kp*y2-rhs];
Print["diffeq2"]
FullSimplify[ diffeq2 /. {alp2->Sqrt[ph/(2*dd)+Sqrt[(ph/(2*dd))^2+kp/dd]],bet2->Sqrt[-ph/(2*dd)+Sqrt[(ph/(2*dd))^2+kp/dd] ]} ]

Print["delta1*ell"]
delta1=Simplify[ell*Integrate[(D[y1,x])^2,{x,c,ell}]/2 /. {alp1->alp1L/ell,bet1->bet1L/ell,c->cL*ell}]
Print["delta2*ell"]
delta2=Simplify[ell*Integrate[(D[y2,x])^2,{x,0,c}]/2 /. {alp2->alp2L/ell,bet2->bet2L/ell,c->cL*ell}]

Print["bc1"]
bc1=y2 /. x->c;
bc1=bc1 /. {alp2->alp2L/ell,bet2->bet2L/ell};
bc1=Simplify[bc1 /. c-> cL*ell]
Print["bc2"]
bc2=y1 /. x->c;
bc2=bc2 /. {alp1->alp1L/ell,bet1->bet1L/ell};
bc2=Simplify[bc2 /. c-> cL*ell]
Print["bc3"]
bc3=y1 /. x->ell;
bc3=bc3 /. {alp1->alp1L/ell,bet1->bet1L/ell};
bc3=Simplify[bc3 /. c-> cL*ell]
Print["bc4"]
bc4=D[y1,x] /. x->ell;
bc4=bc4 /. {alp1->alp1L/ell,bet1->bet1L/ell};
bc4=Simplify[bc4 /. c-> cL*ell]
Print["bc5"]
bc5=ell*D[y2-y1,x] /. x->c;
bc5=bc5 /. {alp2->alp2L/ell,bet2->bet2L/ell,alp1->alp1L/ell,bet1->bet1L/ell};
bc5=Simplify[bc5 /. c-> cL*ell]
Print["bc6"]
bc6=ell^2*Simplify[D[y2-y1,{x,2}] /. x->c];
bc6=bc6 /. {alp2->alp2L/ell,bet2->bet2L/ell,alp1->alp1L/ell,bet1->bet1L/ell};
bc6=Simplify[bc6/. c-> cL*ell]
sol=Solve[{bc1==d0,bc2==d0,bc3== 0,bc4==0,bc5==0,bc6==0},{b3,b4,b7,b8,b9,b10}];
Print["b3"]
bb3=Part[Simplify[b3/d0 /. sol],1]
Print["b4"]
bb4=Part[Simplify[b4/d0 /. sol],1]
Print["b7"]
bb7=Part[Simplify[b7/d0 /. sol],1]
Print["b8"]
bb8=Part[Simplify[b8/d0 /. sol],1]
Print["b9"]
bb9=Part[Simplify[b9/d0 /. sol],1]
Print["b10"]
bb10=Part[Simplify[b10/d0 /. sol],1]
delta=delta1+delta2;

Print["bc7"]
bc7=ell^2*D[y1,{x,2}] /. x->ell;
bc7=Simplify[bc7 /. {alp1->alp1L/ell,bet1->bet1L/ell, c->cL*ell}]
Print["bc8"]
bc8=Simplify[ell^3*D[y2-y1,{x,3}] /. x->c];
bc8=Simplify[bc8 /. {alp2->alp2L/ell,bet2->bet2L/ell,alp1->alp1L/ell,bet1->bet1L/ell}];
bc8=Simplify[bc8 /. c->cL*ell]
Print["eq1"]
eq1=bc7 /. {alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]],alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]]}
Print["eq2"]
eq2=bc8/. {alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]],alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]]}
Save["tempfile2.mth",{eq1,eq2,bb3,bb4,bb7,bb8,bb9,bb10,delta,rhs,kp}]

<<"tempfile2.mth";
phLD=.
kLD=.
phLDG=8.45;
kLDG=4.90^4;
m=9.699
hd=1000
Print["i    c/L   phL^2/2D   ph/2(kD)^(1/2)  L(k/D)^(1/4) delph/2(kD)^(1/2)   y(0)/d0   energyfrac  (ph/LD)^2-kL^4/D    Imag/Real"]
For[iii=1,iii≤ 51,iii++,cLs=0.684+0.0001*(iii-1);
b3s = bb3 /.     {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
b4s = bb4 /.     {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
b7s = bb7 /.     {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
b8s = bb8 /.     {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
b9s = bb9 /.     {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
b10s = bb10 /. {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
delphDk=6*delta/(hd^2*Sqrt[kLD])  /. {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],
                alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
delphDk = delphDk /. {cL->cLs,b3->b3s,b4->b4s,b7->b7s,b8->b8s,b9->b9s,b10->b10s};
delphDk = delphDk /. {cL->cLs};
delp=6*delta/hd  /. {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
eq1c=eq1-delp /. {cL->cLs,b3->b3s,b4->b4s,b7->b7s,b8->b8s,b9->b9s,b10->b10s};
eq1c=eq1c /. cL->cLs;
eq2c=eq2 /. {cL->cLs,b3->b3s,b4->b4s,b7->b7s,b8->b8s,b9->b9s,b10->b10s};
eq2c=eq2c /. cL->cLs;
sol2=FindRoot[{eq1c==0,eq2c==0},{{phLD,phLDG},{kLD,kLDG}}];
pLs=phLD /. sol2; phLDG=pLs;
kLs=kLD /. sol2; kLDG=kLs;
phDk=pLs/Sqrt[kLs];
delphDk=delphDk /. {phLD->pLs,kLD->kLs};
y0=b3s+b4s-rhs/(d0*kp) /. {cL->cLs,phLD->pLs,kLD->kLs};
energyfrac=(1/3)*delphDk^2*hd^2/m;
rad=pLs^2-kLs;
rat1=Im[pLs]/Re[pLs];
rat2=Im[phDk]/Re[phDk];
rat3=Im[kLs]/Re[kLs];
rat4=Im[delphDk]/Re[delphDk];
rat5=Im[y0]/Re[y0];
rat6=Im[energyfrac]/Re[energyfrac];
rat7=Im[rad]/Re[rad];
err=Max[Abs[rat1],Abs[rat2],Abs[rat3],Abs[rat4],Abs[rat5],Abs[rat6],Abs[rat7]];
Print[iii,",  ",cLs,",  ",Re[pLs],",     ",Re[phDk],",     ",Re[kLs]^0.25,",     ",Re[delphDk],",     ",Re[y0],",   ",Re[energyfrac],",   ",Re[rad],",   ",err,";"]]

b3s=b3s /. {cL->cLs,b10->b10s,cpL->cpLs,phLD->pLs,kLD->kLs};
b4s=b4s /. {cL->cLs,b10->b10s,cpL->cpLs,phLD->pLs,kLD->kLs};
b7s=b7s /. {cL->cLs,b10->b10s,cpL->cpLs,phLD->pLs,kLD->kLs};
b8s=b8s /. {cL->cLs,b10->b10s,cpL->cpLs,phLD->pLs,kLD->kLs};
b9s=b9s /. {cL->cLs,b10->b10s,cpL->cpLs,phLD->pLs,kLD->kLs};
y1=b7*Cos[alp1*x]+b8*Sin[alp1*x]+b9*Cos[bet1*x]+b10*Sin[bet1*x]  /. {b1->b1s,b2->b2s,b3->b3s,b4->b4s,b5->b5s,b6->b6s,b7->b7s,b8->b8s,b9->b9s,b10->b10s};
y1 = y1 /. {alp1->alp1L,bet1->bet1L};
y1= y1 /. {alp1L->Sqrt[phLD+Sqrt[phLD^2-kLD]],bet1L->Sqrt[phLD-Sqrt[phLD^2-kLD]],alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]]};
y2=b3*Cos[alp2*x]+b4*Cosh[bet2*x]-rhs/kp  /. {b3->b3s,b4->b4s,b7->b7s,b8->b8s,b9->b9s,b10->b10s};
y2 = y2 /. {alp2->alp2L,bet2->bet2L};
y2=y2 /. {alp2L->Sqrt[phLD+Sqrt[phLD^2+kLD/(m-1)]],bet2L->Sqrt[-phLD+Sqrt[phLD^2+kLD/(m-1)]]};
y1=y1 /. {phLD->pLs,kLD->kLs,cL->cLs,d0->1}
y2=y2 /. {phLD->pLs,kLD->kLs,cL->cLs,d0->1}
Plot[Re[y2],{x,0,Re[cLs]}]
Plot[Re[y1],{x,Re[cLs],1}]