clearclcme=9.109e-31;h=1.054e-34;ee=1.6e-19;V=6.21*ee;zct=10000E=V/zct*(1:zct);ct1=0;n=0;m=1;d=12e-9;a=9e-9;for ct=1:(zct-1)    t=real(tofindzeros(E(ct),n,m,V,d,a))*real(tofindzeros(E(ct+1),n,m,V,d,a));    if t<0        ct1=ct1+1;        x(ct1)=E(ct);    endendx1=x/ee;x1(1:4)for ct=1:4    E=x(ct);    k=(2*E*me/h^2-(m*pi/d)^2).^0.5;    q=(2*V*me/h^2-k.^2).^0.5;    mm=besselj(n,k*a)/besselk(n,q*a);    xx1=linspace(0,a,10001);    yy1=abs(besselj(n,k*xx1)).^2.*xx1;    I1=a/10000*sum(yy1);    xx2=linspace(a,11*a,10001);    yy2=abs(besselk(n,q*xx2)).^2.*xx2;    I2=mm^2*10*a/10000*sum(yy2);    A(ct)=1/(I2+I1).^0.5;    B(ct)=A(ct)*mm;endfor ct=1:4    if ct<=3        ct1=ct;ct2=ct+1;    else        ct1=4;ct2=1;    end    E1=x(ct1);E2=x(ct2);    k1=(2*E1*me/h^2-(m*pi/d)^2).^0.5;    q1=(2*V*me/h^2-k1.^2).^0.5;    k2=(2*E2*me/h^2-(m*pi/d)^2).^0.5;    q2=(2*V*me/h^2-k2.^2).^0.5;        xx1=linspace(0,a,10001);    yy1=conj(A(ct1)*besselj(n,k1*xx1)).*A(ct2).*besselj(n,k2*xx1).*xx1.^2;    I1=a/10000*sum(yy1);    xx2=linspace(a,11*a,10001);    yy2=conj(B(ct1)*besselk(n,q1*xx2)).*B(ct2).*besselk(n,q2*xx2).*xx2.^2;    I2=10*a/10000*sum(yy2);    I=I1+I2;    u(ct)=I;endN=5e24;eer=8.85e-12;w=(1:1000)*2e-5*ee/h;g=1/300e-15;w41=(x(4)-x(1))/h;w31=(x(3)-x(1))/h;w21=(x(2)-x(1))/h;Xi=N*ee^4/eer/h^3*prod(u)./(w41-3*w-i*g)./(w31-2*w-i*g)./(w21-w);hw=h*w/ee;semilogy(hw,abs(Xi))hold on