n=100;
%Zufallszahlen
X=rand(n,2);
%Punkte innerhalb des Kreises
Y1=(X(:,1).^2+X(:,2).^2 <=1);
%Punkte außerhalb des Kreises
Y2=ones(n,1) - Y1;
%Annäherung an pi
Z=4*sum(Y1)/n;
disp(['Pi ist ungefähr ', num2str(Z),'.'])
%plot
hold on
x=0:.01:1;
y=sqrt(1-x.^2);
plot(x,y)
rot=X.*[Y1 Y1];
blau=X.*[Y2 Y2];
plot(rot(:,1),rot(:,2),'r.')
plot(blau(:,1),blau(:,2),'b.')
legend('Einheitskreis','innere Punkte','äußere Punkte')