Question: How can I get the variational value of omega?

Hello:
    I  want to get the variational value of omega when KR[N]>0.2, but  the omega of my procedure remain the same,  how do I change my program to get the variational value of omega when KR[N]>0.2.

Thank you !

My procedure:
restart:
> Digits:=30;
>
> h0:=0.156;
> d:=0.32*h0;
> l:=2;
> h1:=h0-d;
> h2:=h0+d;
> h3:=0.6*h0;
> g:=9.8;
> d1:=1;
> Term:=8;
> Num:=100:
> n:=2:
> l1:=l/n;
>
> for N from 1 to Num do
> lambda:=2*n*Pi/l:
> k0:=evalf(10^(-10)*Pi+2.5*(N-1)*Pi/(Num-1)):
> tau0:=evalf(k0*h0):
> omega:=evalf((g*k0*tanh(k0*h0))^(1/2)):
> E:=g/(omega)^2:
> k1:=abs(fsolve(omega^2=g*k*tanh(k*h1),k)):
> tau1:=evalf(k1*h1):
> k2:=abs(fsolve(omega^2=g*k*tanh(k*h2),k)):
> tau2:=evalf(k2*h2):
> k3:=abs(fsolve(omega^2=g*k*tanh(k*h3),k)):
> tau3:=evalf(k3*h3):
>
>
>
> u0:=evalf(g*tanh(tau0)*(1+2*tau0/(sinh(2*tau0)))/(2*k0));
> u01:=evalf(g*tanh(tau3)*(1+2*tau3/(sinh(2*tau3)))/(2*k3));
> u1:=evalf(g*(1-(tanh(tau0))^2)*(sinh(2*tau0)-2*tau0*cosh(2*tau0))/(4*(2*tau0+sinh(2*tau0))));
> H:=evalf((1+2*tau0/sinh(2*tau0))/(-lambda*k0*d));
> delta00:=evalf((lambda*d*u1/u0+I*k0)*H);
> delta01:=evalf((lambda*d*u1/u0-I*k0)*H);
> delta11:=evalf((lambda*d*u1/u0+I*k0)*H*exp(I*k0*l));
> delta12:=evalf((lambda*d*u1/u0-I*k0)*H*exp(-I*k0*l));
> delta21:=evalf(exp(I*k0*(l+l1)));
> delta22:=evalf(u0*k0*exp(I*k0*(l+l1)));
> delta31:=evalf(exp(-I*k0*(l+l1)));
> delta32:=evalf(-u0*k0*exp(-I*k0*(l+l1)));
> delta41:=evalf(exp(I*k3*(l+l1)));
> delta42:=evalf(u01*k3*exp(I*k3*(l+l1)));
> delta51:=evalf(exp(-I*k3*(l+l1)));
> delta52:=evalf(-u01*k3*exp(-I*k3*(l+l1)));
> delta61:=evalf(exp(I*k3*(l+l1+d1)));
> delta62:=evalf(u01*k3*exp(I*k3*(l+l1+d1)));
> delta71:=evalf(exp(-I*k3*(l+l1+d1)));
> delta72:=evalf(-u01*k3*exp(-I*k3*(l+l1+d1)));
> delta81:=evalf(exp(I*k0*(l+l1+d1)));
> delta82:=evalf(u0*k0*exp(I*k0*(l+l1+d1)));
>
>
> Hi:=(Matrix([[1,1],[delta00,delta01]]))^(-1);
> H1:=(Matrix([[delta21,delta31],[delta22,delta32]]))^(-1);
> H2:=Matrix([[delta41,delta51],[delta42,delta52]]);
> H3:=(Matrix([[delta61,delta71],[delta62,delta72]]))^(-1);
> H4:=Matrix([[delta81],[delta82]]);
>
>
>
>
>
> MM:=evalm(H1.H2.H3.H4):
> R:=evalf(MM[2,1]/MM[1,1]):
> Kr:=abs(evalf(R)):
> KR[N]:=abs(Kr);
>
> end do:
>
> L:=seq(KR[N],N=1..Num):
>
> for r  in L do
> if r>0.2 then
> print(omega);
> end if;
> end do;

Please Wait...