N=100000;
Ncenter=0; %counter
Nedge=0; %counter
for i=1:N
    x=rand;
    y=rand;
    D2toCenter=(x-0.5)^2+(y-0.5)^2;
    D2toEdge=min([x^2 y^2 (1-x)^2 (1-y)^2]);
    if (D2toCenter<D2toEdge)
        Ncenter=Ncenter+1;
    else
        Nedge=Nedge+1;
    end
end
area=Ncenter/N