
% Probability density function: natural measue
% Ricardo Carretero Nov'22

%To see evolution:
%clear;for Nit=2.^[1:17];pdf;end

fs=15;
lw=2;

f = inline('4*x.*(1-x)','x');
if(exist('Nit')==0) Nit=1e4;end
Nbox=125;
Ntrans=1000;
x=zeros(Nit,1);

figure(1);clf

xo=0.01;
for ii=1:Ntrans % iterate transients
 xn=f(xo);
 xo=xn;
end

for ii=1:Nit % keep Nit iterates
 xn=f(xo);
 xo=xn;
 x(ii)=xn;
end

h=histogram(x,Nbox,'Normalization','pdf');
set(gca,'FontSize',[fs]);

xx=linspace(0,1,500);
hold on
plot(xx,1./(pi*sqrt(xx-xx.^2)),'LineWidth',lw);
hold off

axis([0 1 0 5]);
xlabel('x')
ylabel('p(x)');
title(['Nit=',num2str(Nit)])


