b-s公式倒是明白，但是没用来计算过 credit risk，楼主的程序我给帮忙改了改，不知道是不是楼主的意思了
function[Va,SigmaVa]=CalV(E,SigmaE,D,r,T)
%a=zeros(8,8);
D=[
    1096853709
    1422513664
    1080412570
    1099330560
    821720
    676976
    196840
    37877619
    ];
a=D;
%b=zeros(8,8);
E=[
    5482334
    112522064
    57679892
    4537800
    751520
    1882103.375
    285241.125
    3833159
    ];
b=E;
%c=zeros(8,8);
SigmaE=[
    0.401835787
    0.543863954
    0.609681884
    0.741317746
    0.917986928
    0.931843334
    1.004589468
    1.02883818
    ];
c=SigmaE;
%计算Va'SigmaVa
key=0;
r=0.0225;
T=1;
PI=4*atan(1);
Va=b;%va以E为迭代的初值
newVa=Va;
SigmaVa=c; %SigmaVa以SigmaE为迭代的初值
for k=1:10  %迭代SigmaVa
    for j=1:10 %迭代Va
        oldVa=newVa;
        % Va=(E+D*exp(-r,T)*Nd2)/Ndl;
        % SigmaVa=E*SigmaE/(E+D*exp(-r,T)*Nd2);
        %t1=log(oldVa/a);
        d1=(log(oldVa./a)+(r+(SigmaVa.^2)/2)*T)./(SigmaVa.*(T^(1/2)));
        d2=d1-(SigmaVa*(T^(1/2)));
        Nd1=normcdf(d1,0,1);
        Nd2=normcdf(d2,0,1);
        fV=oldVa.*Nd1-a.*exp(-r*T).*Nd2-b;
        temp111=exp(-d1.^2/2)./((2*PI*T)^0.5.*SigmaVa)
        dfV=Nd1+temp111-a.*exp(-r*T).*exp(-d2.^2/2)./((2*PI*T)^0.5.*oldVa.*SigmaVa);
        if dfV==0
            dfV
            pause;
        end
        dV=-fV./dfV;
        newVa=oldVa+dV;   %前后两次的值认为非常接近，则迭代成功
        newVa
        if newVa==0
            oldVa
            newVa
            pause;
        end
        if(dV/newVa)<1e-7 & (dV/newVa)>-1e-7
            break;
        end
    end
    Va=newVa;
    d1=(log(Va./a)+(r+(SigmaVa.^2)/2)*T)./(SigmaVa.*(T^(1/2)));
    d2=d1-(SigmaVa*(T^(1/2)));
    Nd1=normcdf(d1,0,1);
    Nd2=normcdf(d2,0,1);  %方程右边与左边的比值
    diff=(Va.*Nd1-a.*exp(-r*T).*Nd2)./b;
    if diff<1-1e-5 | diff>1+1e-5
        assert('迭代出现错误'); %迭代出现错误
    else
        key=key+1;
    end
    newSigmaVa=b.*c./(b+a.*exp(-r*T).*Nd2);
    if((SigmaVa-newSigmaVa)./newSigmaVa)<1e-5 & ((SigmaVa-newSigmaVa)./newSigmaVa)>(-1e-5)
        break;
        key=key+1;
    else
        SigmaVa=newSigmaVa;
    end
end
if key==2
    %else
    % Va=一l：
    % SigmaVa=-1：
end

烦请楼主多看看基础的matlab知识，还有就是那两个for 循环，其实没有用到，

我把运行的结果给贴在这了，烦请楼主看看正确与否

ans =

  1.0e+009 *

                 NaN
   1.502597151647248
   1.113130803312277
   1.079132735700452
   0.001534193015389
   0.002534900899513
   0.000471079307586
   0.039933626244482

哇，大师，你太牛了。送你500论坛币。请你在这个链接里面回复任何东东，我就把悬赏币给您。https://bbs.pinggu.org/thread-1285211-1-1.html
不过还有个小问题，程序倒是能算出来Va，可是SigmaVa为什么就没有算出来呢？

太兴奋了，哈哈，谢谢大师啊，送你500论坛币吧。请你在这个链接里面回复任何东东，我就把悬赏币给您。https://bbs.pinggu.org/thread-1285211-1-1.html

加菲猫大师，那您改过的这段代码里面，两个for循环都用到了吧？（我不清楚你为什么说我之前的程序没有用到for循环。）

果需要避免NaN，则需要分析一下NaN产生的原因，再看看如何在计算之初就避免NaN出现
有时候出现了NaN，也可以在后续计算中将NaN数据排除在外，比如, MATLAB很多函数 sum、max、min、mean 等都提供'omitnan'的选项，利用这些函数在后续计算里忽略调NaN