wd1e1.m

用二分法

clear;
global wd0 wd1 ws1 ws2  mw0 mb0 ab aw U0 h ;
U0=0.225;c=2.997925*1e+8;
me=0.5109989*1e+6/(c^2);
mw0=0.067*me;mb0=0.092*me;
wd0=85e-10;
ws1=80e-10;
ws2=30e-10;
tev=(1e+19)/1.602177;
h=tev*1.05457266e-34;

aw=0.642;
ab=mw0*aw/mb0;
pts=30;dw=1e-11;
E0=0.02;
%E0=0.047450376;
%E0=0.17804851;
de=0.001;ps=(0.035094804-E0)/0.001;E1=[];E2=[];
b3=[];x=[];q1=[];q2=[];
for j=1:pts
    wd1=(j-1)*1e-10;
    aenerg1=[];
           for k=1:ps
        e1=E0+(k-1)*de;
        e2=e1+de;
         b1=eqs144(e1);
         b2=eqs144(e2);
        if b2*b1<0
            b3=fzero('eqs144',[e1,e2]);
            aenerg1=[aenerg1;b3];
           
        end
        if length(aenerg1)>1
            break
        end
    end
    if length(aenerg1)<1
        aenerg1=[nan;nan];
    elseif length(aenerg1)<2
        aenerg1=[aenerg1;nan];
    else
        aenerg1=aenerg1;
    end
     E1=[E1;aenerg1(1,1)];
     E2=[E2;aenerg1(2,1)];
      x=[x;wd1];
end
save E1wd1.dat E1 -ascii;
plot(x,E1,'k');hold on;
xlabel('wd1');hold on;
ylabel('Energy E1(ev)');hold on;
%save xwd1.dat x -ascii;
return