Conditional Mean simulation

Question

The joint p.d.f of $f(x,y)=2 $ if $0<x<y<1$ and zero else where, I get the conditional expectation $\mathbb{E}(Y|X=x)=\dfrac{x+1}{2}, ~0<x<1$, I verified this with my Matlab code

 for i=1:500000
 x(i)=rand;
 y(i)=x(i)  + (1-x(i))*rand(1,1);
end
plot(x,y)
hold on;  
my=0.5*(x+1);
plot(x,my)
hold on
p=polyfit(x,y,1);
pv=polyval(p,x);
plot(x,pv)

I get the desired result. My question is (i) Is my logic correct when I say that these mathematical answer and Matlab's polyfit match when there INDEED is a linear relationship. (ii) How do I simulate on Matlab an example where $f(x,y)=kxy$