I would like to make surface plot of my probability of failure curve?

Below is the main file and my Function Code an the bottom, i wanna plot ptobability of failure with changing the muv (velocity) and mub (pier width) in the same and plot in surface plot. I attached the function code (scourdepth1) and main code file (POFPOF)

clear all;
mue=0.82;                       %mean of unbiased factor  of study case 
sige=0.189;                     %SD of unbiased factor of study case  
muk1=0.9;                       %lower limit of K1 of study case 
sigk1=1.1;                      %upper limit K1  of study case 
muk2=1;                         %lower limit of K of study case 
sigk2=1;                        %upper limit of K2 of study case 
muk3=1.1;                       %lower limit  of K3 of study case 
sigk3=1.3;                      %upper limit of K3 of study case 
mud50=0.01;                     %Mean of D50
sigd50=0.0195;                  %SD of D50
muv=2.48;                      %mean value of stream velocity 
sigv=0.28;                      %SD of stream velocity COV=0.28 due to Time Variant
muy=4.042;                      %mean value of stream depth of study case 
sigy=0.4042;                    %SD of stream depth of study case COV=0.1 due to Time Variant
%mub=1.067;                      %mean of pier width of study case 
sigb=0.0534;                    %SD of pier width of study case 
i=1;
for mub=0.5:0.1:5
    [Ds] = scourdepth1(muv,sigv,muy,sigy,mub,sigb,muk1,sigk1,muk2,sigk2,muk3,sigk3,mue,sige,mud50,sigd50);
    count=nnz((3-Ds(:))<0);
    total=size(Ds,2);
    b1(i)=mub;
    p(i)=count/total;
    i=i+1;
end
plot(b1,p,'red')
title('Probability of Failure Vs. Flow Velocity');
xlabel('Flow Velocity (m/sec)');
ylabel('Probability of Failure');
grid on

function [Ds] = scourdepth1(muv,sigv,muy,sigy,mub,sigb,muk1,sigk1,muk2,sigk2,muk3,sigk3,mue,sige,mud50,sigd50)
      
for i=1:10000
e(i)=normrnd(mue,sige);
k1(i)=unifrnd(muk1,sigk1);                     %correction factor for pier shape
k2(i)=unifrnd(muk2,sigk2);                     %correction factor of angle of attack
k3(i)=unifrnd(muk3,sigk3);                     %correction factor of bed configuration
D50(i)=lognrnd(mud50,sigd50);                  %Median Grain Size
v(i)=normrnd(muv,sigv);                        %flow velocity
y(i)=lognrnd(muy,sigy);                        %flow stream depth
b(i)=lognrnd(mub,sigb);                        %effective pier width
k4(i)=0.35*((b(i)/D50(i))^0.19);               %Correction Factor for soil properties grain size factor  
f(i)=v(i)./(((9.81.*y(i))^0.5));               %Froude number
Ds(i)=e(i).*2.*y(i).*k1(i).*k2(i).*k3(i).*k4(i).*((b(i)./y(i))^0.65).*(f(i)^0.43);    %maximum Predicted scour depth
end
end

Answer 1

[Torsten]

clear all;
mue=0.82;                       %mean of unbiased factor  of study case 
sige=0.189;                     %SD of unbiased factor of study case  
muk1=0.9;                       %lower limit of K1 of study case 
sigk1=1.1;                      %upper limit K1  of study case 
muk2=1;                         %lower limit of K of study case 
sigk2=1;                        %upper limit of K2 of study case 
muk3=1.1;                       %lower limit  of K3 of study case 
sigk3=1.3;                      %upper limit of K3 of study case 
mud50=0.01;                     %Mean of D50
sigd50=0.0195;                  %SD of D50
%muv=2.48;                      %mean value of stream velocity 
sigv=0.28;                      %SD of stream velocity COV=0.28 due to Time Variant
muy=4.042;                      %mean value of stream depth of study case 
sigy=0.4042;                    %SD of stream depth of study case COV=0.1 due to Time Variant
%mub=1.067;                      %mean of pier width of study case 
sigb=0.0534;                    %SD of pier width of study case 
j=1;
for muv = 0.5:0.5:5
    i=1;
for mub=0.5:0.5:5
    [Ds] = scourdepth1(muv,sigv,muy,sigy,mub,sigb,muk1,sigk1,muk2,sigk2,muk3,sigk3,mue,sige,mud50,sigd50);
    count=nnz((3-Ds(:))<0);
    total=size(Ds,2);
    b1(i,j)=mub;
    b2(i,j)=muv;
    p(i,j)=count/total;
    i=i+1;
end
j=j+1;
end
surf(b1,b2,p)

function [Ds] = scourdepth1(muv,sigv,muy,sigy,mub,sigb,muk1,sigk1,muk2,sigk2,muk3,sigk3,mue,sige,mud50,sigd50)
      
for i=1:10000
e(i)=normrnd(mue,sige);
k1(i)=unifrnd(muk1,sigk1);                     %correction factor for pier shape
k2(i)=unifrnd(muk2,sigk2);                     %correction factor of angle of attack
k3(i)=unifrnd(muk3,sigk3);                     %correction factor of bed configuration
D50(i)=lognrnd(mud50,sigd50);                  %Median Grain Size
v(i)=normrnd(muv,sigv);                        %flow velocity
y(i)=lognrnd(muy,sigy);                        %flow stream depth
b(i)=lognrnd(mub,sigb);                        %effective pier width
k4(i)=0.35*((b(i)/D50(i))^0.19);               %Correction Factor for soil properties grain size factor  
f(i)=v(i)./(((9.81.*y(i))^0.5));               %Froude number
Ds(i)=e(i).*2.*y(i).*k1(i).*k2(i).*k3(i).*k4(i).*((b(i)./y(i))^0.65).*(f(i)^0.43);    %maximum Predicted scour depth
end
end

[Khaled aL Qattawi]

thank you, that is perfect.