The reason is because of the function RHO3. The input "x" is a vector of "x"-coordinates. Thus, in the code, "y" is also assigned to 1 only.
The output "y" ( the values of "d") should be a vector, as well.
Try the following code:
function y=rho3(x)
% radial density function 3
y =ones(size(x));
idx = find((x < 0.5) & (x >= 0.0));
y(idx) = 2
