clear all
close all

% dati del problema
a = input('inserisci il radicando ');
n = input('inserisci l''indice di radice ');
% if mod(n,2)==0 && a<0, disp('impossible'); return; end

f=@(x) x.^n-a;
figure; fplot(f); hold on; yline(0); grid on; hold off
% 
% % iterato iniziale
x0 = input('inserisci approssimazione iniziale ');

% criteri di arresto
% tol = input('inserisci la tolleranza per il criterio di arresto ');
tol = 1e-6;
kmax = 100;
rex = a^(1/n);

[x,k,xvec] = newton_radici(n,a,x0,kmax,tol);
if isreal(x)
    fprintf('\nRadice calcolata: %4.16f\n',x);
else
    fprintf('\nRadice calcolata: %4.16f+%4.16f i\n',real(x), imag(x));
end
fprintf('Numero di iterazioni=%d \n', k);
fprintf('Errore=%d \n', abs(x-rex));

if a>=0
    figure; plot(1:k,xvec,'*r',1:k,rex*ones(size(xvec)),'b-'); 
    xlabel('Iterazioni');
    legend('Approssimazioni','valore esatto')
    title('Approssimazioni');
end


figure; plot(abs(diff(xvec)),'-*'); title('|x_{k+1}-x_k|');

% -------------------------------------------------------------------------
function [x,k,xvec] = newton_radici(n,a,x0,kmax,tol)
% newton_radici: calcola la radice n-esima di un numero
% INPUT: n = indice del radicale
%        a = numero del quale si vuole determinare la radice quadrata 
%        b = approssimazione iniziale
%        kmax = massimo numero di iterazioni
%        tol = tolleranza per il criterio di arresto sulla differenza tra
%        due approssimazioni successive
% OUTPUT: x = approssimazione calcolata
%         k = numero di iterazioni
%         bvec = vettore delle approssimazioni successive



% Inizializza 
xvec = zeros(kmax+1,1);
xvec(1) = x0;

fprintf('k=1, x=%4.8f \n', xvec(1));
if x0==0
    disp('Il metodo non può iniziare'); 
    x=x0; k=1; xvec=x0;
    return; 
end

k=1;
continua = true; 
while continua 
    
    % aggiorna la stima
    xvec(k+1) = xvec(k)-(xvec(k)^n-a)/(n*xvec(k)^(n-1));
    
    % visualizza le stime calcolate
    if isreal(x0)
        fprintf('k=%d, x=%4.16f, err=%4.16f \n', k+1, xvec(k+1),abs(xvec(k+1)-a^(1/n)));
    else
        fprintf('k=%d, x=%f+%f i, err=%4.16f \n',k+1, real(xvec(k+1)),imag(xvec(k+1)),abs(xvec(k+1)-a^(1/n)));
    end

    % controlla condizioni di arresto 
    diff = abs(xvec(k+1)-xvec(k));
    if xvec(k+1)==0
        disp('Il metodo non può continuare');
        continua = false;
    else
        continua = k<kmax & diff >tol;
    end

    k = k+1;

end
xvec = xvec(1:k);
x = xvec(end);
end




    