我的公式和公式说明都在里面了。自己试着编了下matlab，但是结果不对啊，现在我调整了下语句，变成fortran语句了，尽管有些可能不大地道，大家帮忙看下这样对不对啊。难啊！！苦啊！！累啊！！
do Ef=-1.0,1.0,0.021
T=0; 
tt=0;%储存变量，j<tt，D(j,j)<Ef,D(j+1,j+1)>Ef
%先判断费米能在特征值里处在什么位置，tt
do j=1,N-1
    if (D(j,j)<Ef)&&(D(j+1,j+1)>Ef) then
         tt=j;
    End if
end do
%费米能以上和以下的特征值就分开了，按顺序参与运算，求和
 do m=1,tt
        S=0;
      do k=tt+1:N
        Sum1=V(:,m)'*Vx*V(:,k)*V(:,k)'*Vy*V(:,m);
        Sum2=conj(Sum1);
        Sum=-1i*(Sum1-Sum2)/(D(m,m)-D(k,k))^2;
        S=S+Sum;
      end do
      T=S+T;
 end do
 ！下面是matlab里的画图语句
 plot(Ef,T)
 hold on
 
end do
保证每一个Ef对应一个T，最后画出Ef--T图形。