function y=f(n,x,r) 
  y(1,1)=x(1); y(1,2)=x(2)
  for i=2:n
    y(i,1)=y(i-1,2); y(i,2)=r*y(i-1,2)*(1-y(i-1,1))
  end
endfunction

ini=["2.100","00200"];
[ok,r,n]=getvalue("Dare il valore di r e n",["r","n"],..
  list("col",1,"col",1),ini);

xset("window",0);         
xbasc(); 
labels=["xmin";"xmax";"ymin";"ymax"];
descr="Estremi rettangolo coordinate";
xmin=0;xmax=1;ymin=0;ymax=1;
plot2d(0,0,strf="015",rect=[xmin ymin xmax ymax]);
xtitle("PIANO DELLE FASI");
//x_message('assegnare con un click la condizione iniziale');
[c_i,c_x,c_y]=xclick();  
xin=[c_x;c_y];    // assegno la condizione iniziale cliccandovi sopra con il mouse 
risp=1;    // introduco una variabile di scelta
 
while %t
  if ((risp==6)|(risp==0)) then break;
  end
  
  while risp==1
    y=f(n,xin,r);
    plot2d(y(:,1),y(:,2),style=0,strf="115",rect=[xmin ymin xmax ymax],leg="ORBITA");
    xtitle("PIANO DELLE FASI");
    xin=y($,:);      
    risp=x_choose(["Continua","Cancella e continua",..
       "Cancella e cambia condizione iniziale","Aggiunge nuova traiettoria","Nuove coordinate",..
       "Fine"],"Cosa fare?");    // l'utente sceglie tra le opzioni
    
    // l'opzione 2 consiste nel graficare la soluzione solo per un intervallo
    // di tempo di ampiezza n cancellando il grafico precedente
    if risp==2 then 
     xbasc();
     plot2d(0,0,strf="115",rect=[xmin ymin xmax ymax],leg="ORBITA")
     xtitle("PIANO DELLE FASI");
     risp=1; 
    end
    
    if risp==3 then 
     xbasc();
     plot2d(0,0,strf="115",rect=[xmin ymin xmax ymax],leg="ORBITA")
     xtitle("PIANO DELLE FASI");
     [c_i,c_x,c_y]=xclick();  
     xin=[c_x;c_y];    // assegno la condizione iniziale cliccandovi sopra con il mouse 
     risp=1; 
    end

    if risp==4 then 
     [c_i,c_x,c_y]=xclick();  
     xin=[c_x;c_y];    // assegno la condizione iniziale cliccandovi sopra con il mouse 
     risp=1; 
    end

    // l'opzione 5 permette di cambiare le coordinate del piano delle fasi per meglio visualizzare il grafico
    if risp==5 then 
     xbasc();
     ini=[string(xmin);string(xmax);string(ymin);string(ymax)];
     [ok,xmin,xmax,ymin,ymax]=getvalue(descr,labels,..
       list("col",1,"col",1,"col",1,"col",1),ini);
     plot2d(0,0,strf="115",rect=[xmin ymin xmax ymax],leg="ORBITA");
     xtitle("PIANO DELLE FASI");
     xin=y(:,1);
     risp=1; 
    end
    
  end
end
