clear
clf
c=input(' Give value of cos alpha: ');
s=sqrt(1-c^2);
nmax=input(' Give number of applications of mapping: ');
npep=input(' Input peppering density parameter i.e. 2/delx (e.g. 5): ');
saxis=1;
delx=2/npep; % might want to change for some cos-alpha values
R2=20;
x=-1:delx:1; % Create grid of initial points; might want to change (see above)
y=x;
[xx,yy]=meshgrid(x,y);
hold on
for i=1:nmax
   temp=xx*c-(yy-xx.*xx)*s; %Do mapping
   xxnew = min(abs(temp),R2).*sign(temp);
   yynew=xx*s+(yy-xx.*xx)*c;
   plot(xxnew,yynew,'.','MarkerSize',4)
   xx=xxnew;
   yy=yynew;
end
while saxis ~= 0
 axis([-saxis saxis -saxis saxis]);
 saxis=input(' enter axis-scaling, usually 1, 0 to stop: ');
end
axis square