Hi, I've written a MATLAB code for this problem, but I'm getting terribly garbage answers (the correct answer is 6). What am I doing wrong?

(The problem is attached)

clear all;
clc;
x(1) = rand();
y(1)=rand();
n= 10^6;
s=x(1) + y(1);
dev2=0;
for i = 1:n-1
rho(i)=rand();
%t(i)= - (reallog (rho(i)/2))/4;
% p(i)= - (reallog (rho(i)/0.4))/0.8;
if rho(i)0.5
t(i) = -(reallog(2-2*rho(i)))/4;
p(i) = -(reallog(2-2*rho(i)))/0.8;

end
fxy(i) = (1/(pi*(1-0.5^2)))*exp(-((x(i)-5)^2 + ((y(i)-1)^2)/0.5^2 - 2*(x(i)-5)*(y(i)-1))/(2*(1-0.5^2)));
ftp(i) = (1/(pi*(1-0.5^2)))*exp(-((t(i)-5)^2 + ((p(i)-1)^2)/0.5^2 - 2*(t(i)-5)*(p(i)-1))/(2*(1-0.5^2)));
gu(i) = 2*exp(-4*abs(t(i)));
gx(i) = 2*exp(-4*abs(x(i)));
hv(i) = 0.4*exp(-0.8*abs(p(i)));
hy(i) = 0.4*exp(-0.8*abs(y(i)));
R(i) = (ftp(i)*hy(i)*gx(i))/(fxy(i)*hv(i)*gu(i));
m(i) = min(1,R(i));
if rho(i)