delta=0.75;
for i=2:n
APHNEW=expm([(i-1).*H,(i-1).*(2*(Q'*Q));(i-1).*(R+M),(i-1).*(-(H'))]);
APHNEW11=APHNEW(1:2,1:2);
APHNEW21=APHNEW(3:4,1:2);
APHNEW12=APHNEW(1:2,3:4);
APHNEW22=APHNEW(3:4,3:4);
BPHNEW22=inv(APHNEW22);
PSIi=BPHNEW22*APHNEW21;
SPSI1=SPSI1+PSIi;
PHIi=beta*(log(det(APHNEW22))+(i-1)*trace(H'))/2;
S0i(i)=exp(-((rbar+mubar)*(i-1)+PHIi));
S0i2(i)=(-((rbar+mubar)*(i-1)+PHIi));
Yn0=Yn0+S0i2(i);
fun0=@(Gamma) exp(trace(1i*(THETA1*inv(eye(2)-2i*SIGMA*Gamma)*SIGMA*Gamma)))/((det(eye(2)-2i*SIGMA*Gamma))^(beta/2));
fun1 = @(Gamma1) exp((1i*Gamma1+delta)*Yn0)*S0i(i)*fun0(1i*PSIi-(Gamma1-1i*delta)*SPSI);
fun3 = @(Gamma1) exp((1i*Gamma1+delta)*Yn0)*(S0i(i)*fun0(1i*PSIi-(Gamma1-1i*delta)*SPSI)-K*fun0(-(Gamma1-1i*delta)*SPSI));
end