eta=-1; k=0.01; ro=i*eta/2/(1-i*eta/2); to=1/(1-i*eta/2); for in=0:99 k=2*pi/100*in; p=exp(i*k); R=ro; T=to; dat=[]; for j=1:100 r=ro+(rand(1)+i*rand(1))/100.*0; t=to+(rand(1)+i*rand(1))/100.*0; Rn=R+T^2*r*p^2/(1-p^2*r*R); Tn=T*t*p/(1-p^2*r*R); R=Rn; T=Tn; abs(R)^2+abs(T)^2; dat=[dat ; T]; end hold off plot(dat,'b-') axis([-1 1 -1 1]) hold on plot(real(dat(1)),imag(dat(1)),'r*') in pause % if (in<10) % eval(['print -dgif8 out0' num2str(in) '.gif']); % else % eval(['print -dgif8 out' num2str(in) '.gif']); % end end