function nes = nestabilni_metoda_REKURZE(usesingleprec)

% Demonstrace nestabilni rekurze Phi(n) = Phi^n, kde Phi = (-1+sqrt(5))/2
%   a rekurzni vzorec je Phi(n+1) = Phi(n-1) - Phi(n)
% Nestabilita vznika proto, ze 2.koren rovnice -(1+sqrt(5))/2 
%   je v absolutni hodnote vetsi nez dany koren
% Nastaveni parametru usesingleprec=1 emuluje jednoduchou presnost
%

%% To which power we are going
niter = 50;

%% Arrays for saving of the results
Phin_multi = zeros(niter,1);
Phin_recur = zeros(niter,1);
err_rel    = zeros(niter,1);

disp(['  '])
if (usesingleprec==1)
    disp(['Emulating single precision'])
    Phi    = single(0.5*(sqrt(5)-1));
    fisouc = single(1);
else
    disp(['Using default Matlab precision'])
    Phi    = 0.5 * (sqrt(5) - 1);
    fisouc = 1;
end
disp(['  '])

%% Initialize
finew  = Phi;
fipr   = Phi;

%whos fisouc 
%whos finew 
%whos fipr 

Phin_multi(1) = fipr;
Phin_recur(1) = finew;
err_rel(1)    = (finew-fipr)/fipr;

%% Do the calcualtion
for n = 2:niter
    % Calculating in working precision
    fistar  = fisouc;
    fisouc  = finew;
    fipr    = fipr * Phi;
    finew   = fistar - fisouc;
    
    % Saving to default precision arrays
    Phin_recur(n) = finew;
    Phin_multi(n) = fipr;
    err_rel(n)    = (finew - fipr)/fipr;
end

%whos  Phin_recur
%whos  Phin_multi
%whos  err_rel    

%% Plot it
figure(1+usesingleprec)
power = (1:niter);
pm = plot(power,Phin_multi,'rx','linewidth',2);
hold on;
pr = plot(power,Phin_recur,'bo-','linewidth',2);
pe = plot(power,err_rel,'g--','linewidth',2);
hold off;
axis([1 niter -1 1])

leg=legend([pm pr pe],...
           {'multiplication';'subtraction';'rel. diff. (error)'},...
           'location','sw');
xlabel('EXPONENT n');
ylabel('CALCULATED \Phi^n');
if (usesingleprec==1) 
    title('SINGLE PRECISION');
else
    title('DOUBLE PRECISION');
end;


%% Print it
%disp([' '])
%format long
%disp(['   Phi^n multip.       Phi^n recur.        rel.err.'])
%[Phin_multi Phin_recur err_rel]