我用二分查找法解一个带有变限二重积分的二元方程组，（未知数x1，x2），蓝色字体的计算式是正确的，但是我不知道怎么把他们放在循环语句里，所以导致运行不出来，哪位大神帮忙看一下哪里有错误，或者怎样做可能不易出现错误？
方程1：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=0
方程2：Phi2=30*normcdf(x2,mu2,sigma2)+90*P2+90*P1+70*P4+90*P3+70*P5-55=0



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