function y=reliabi(s0001,s0002)
fprintf('\n')
fprintf('It is running reliabi,please waiting......'); 
%定義二組隨機變量  
% define 2 randon varibles for colounm(s01,s02)
% corresponding equivalent randon varible colounm(ys01,ys02)
syms x1     x2     x3    x4    x5   x6   x7   H    c     a   b   fa
syms ux1    ux2    ux3   ux4   ux5  ux6  ux7  uH  uc     ua  ub  ufa
syms fx1    fx2    fx3   fx4   fx5  fx6  fx7  fH  fc     fa  fb  ffa
s01=[x1     x2     x3    x4    x5   x6   x7   H    c     a   b   fa];
s001=[ux1    ux2    ux3   ux4   ux5  ux6  ux7  uH  uc   ua  ub  ufa];
s002=[fx1    fx2    fx3   fx4   fx5  fx6  fx7  fH  fc   fa  fb  ffa];
%considering element load and node load ,randon varible is x01,corresponding average value and 
%方差 are s001,s002.
ini_Ks=init_ka   %抗滑移、抗傾覆和承載力參數
x10=s0001;
y10=(x10-s0001)./s0002;y10=1e-5*round(1e5*y10);r1m=sqrt(y10*y10')
yequival=G(y10,x10,s0001,s0002);
y20=yequival{1,1};r3m=sqrt(y20*y20')
x20=s0002.*y20+s0001;x20=1e-5*round(1e5*x20);
m=1; 
while ((abs(r1m-r3m)>=1e-2) & (m<15))
    x10=x20;r1m=r3m;
    y10=(x10-s0001)./s0002;y10=1e-5*round(1e5*y10);r1=sqrt(y10*y10');   
    yequival=G(y10,x10,s0001,s0002);
    y20=yequival{1,1}; r3m=sqrt(y20*y20');
    if r3m>8  & m>8
        r1m=r3m;
    end
    fprintf('\n');
    fprintf('r1m=%f   r3m=%f    ',r1m,r3m);
    x20=s0002.*y20+s0001;x20=1e-5*round(1e5*x20);m=m+1;
end

beta=r1m;pfm=failprob(beta);
AmZm=yequival{1,2};

beta  

y01=(s01-s001)./s002;
n1=length(s0001);
dy01=ones(1,n1)*jacobian(y01,s01);
dux01=ones(1,n1)*jacobian(y01,s001);
dfx01=ones(1,n1)*jacobian(y01,s002);

dy01=subs(dy01,s01,x20);dy01=subs(dy01,s001,s0001);dy01=subs(dy01,s002,s0002);
dux01=subs(dux01,s01,x20);dux01=subs(dux01,s001,s0001);dux01=subs(dux01,s002,s0002);
dfx01=subs(dfx01,s01,x20);dfx01=subs(dfx01,s001,s0001);dfx01=subs(dfx01,s002,s0002);

duy01=-AmZm'.*dy01;
dux01=-AmZm'.*dux01;
dfx01=-AmZm'.*dfx01;

duy=numeric(duy01/beta);
dux=numeric(dux01/beta);
dfx=numeric(dfx01/beta);
 
y={duy dux dfx beta};
beta
duy'
dux'
dfx'