扩展到二维以后，如果变量之间互相独立，normcdf可以用，看matlab帮助，normcdf允许输入向量或矩阵。
按第一种方法求解，如果用二重积分，只要能写在一个函数中，应该也是可以的。

你只需要把syms变量抛开，这是你最主要的障碍。在自定义函数中，你将自定义函数写成如下形式：

r=-40*mvncdf([x(1) x(2)],mu,sigma)+45*(1-mvncdf([x(1)+150    x(2)],mu,sigma)+…..

我没有注意到有mvncdf函数，有了它，连独立性假设都不要了。
只需按我前面所提供的例子写自定义函数就行了，无论采用哪种方法。注意自定义函数返回一个标量数值就行。

对于5和6，可以考虑变量替换，如x3=x1+x2，这样的话就可以正常使用mvncdf函数了，但它又带来了一个新问题，新的多元变量为[x1 x3]，此时协方差矩阵要重新给出。积分是个十分费力的事。

我用二分查找法解x1和x2，涉及到积分的式子我试过，是正确的，但是放在循环语句里，就运行错误，麻烦您帮我看一下是哪里出现了错误，谢谢！
ZERO=10^-4;
ZERO1=10^-8;
a=150;
mu1=400;sigma1=112;
mu2=300;sigma2=105;
x2_min=0;x2_max=800;
Error=1;
     while(Error&&(x2_max-x2_min)>ZERO)
     x2=(x2_min+x2_max)/2;
     x1_min=0;x1_max=800;
     Error1=1;
        while(Error1&&(x1_max-x1_min)>ZERO1)
        x1=(x1_min+x1_max)/2;
        syms x y
        P1=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x1+x2+a-x,+inf),x,x1,x1+a));
        Phi1=-40*normcdf(x1,mu1,sigma1)+20*( normcdf(x1+a,mu1,sigma1)-normcdf(x1,mu1,sigma1))+45*(1-normcdf(x1+a,mu1,sigma1))+20*P1;  
            if Phi1<-ZERO1
                 x1_max=x1;
                 elseif Phi1>ZERO1
                 x1_min=x1;
                 else
                 Error1=0;
            end
        end
        x1
        syms x y
        P2=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x2,+inf),x,x1+a,+inf));
        P3=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x2+a,+inf),x,0,x1));
        P4=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x2,x1+x2+a-x),x,x1,x1+a));
        P5=vpa(int(int(mvnpdf([x y],[mu1 mu2],[sigma1^2 0;0 sigma2^2]),y,x2,x2+a),x,0,x1));
        Phi2=30*normcdf(x2,mu2,sigma2)+90*P2+90*P1+70*P4+90*P3+70*P5-55;
        if Phi2<-ZERO
        x2_max=x2;
        elseif Phi2>ZERO
        x2_min=x2;
        else
        Error=0;
        end
     end
x1
x2