% NMA_Newton vice iteraci, overeni, vypocet norem
syms x y                  % deklarace symbolickych promennych

f = -1*(x-1)^2 - y + 1;   % zadane rovnice f=0, g=0
g = x^3 - y ;

X0 = [0.5 ; 0.5];         % zvolene X0

F = [f; g];               % vektor F
J=jacobian(F,[x y]);      % matice J


N  = 3;

% budouci vysledky
Xn  = zeros(2,N+1);  % vynulovana matice 2 radky, N+1 sloupcu
rez = zeros(1,N+1);  % vektor N+1 rezidui (norem F(Xn))

Xn(:,1) = X0;


for k = 1 : N        % vypocet X1, ..., XN

    FXn = subs(F,{x y},Xn(:,k)');    % dosazeni Xn
    JXn = subs(J,{x y},Xn(:,k)');

    FXn = double(FXn);               % prevedeni symbolickych vyrazu
    JXn = double(JXn);               % do ciselnych

    % overeni vhodnosti aktualniho Xn
    if ~(any(isreal(FXn)))
         error('nevhodny bod (F)')
    end

    if ~(any(isreal(JXn)))
        error('nevhodny bod (J)')
    end

    if ~(any(isfinite(FXn)))
        error('nevhodny bod (F)')
    end

    if ~(any(isfinite(JXn)))
        error('nevhodny bod (J)')
    end

    if abs(det(JXn)) < eps
        error('nevhodny bod (J)')
    end

    rez(k) = norm(FXn, 2);           % euklidovska norma F(Xn)

    DeltaN = JXn\(-FXn)              % vypocet "Delta", reseni lin. rovnic
    Xn(:,k+1) = Xn(:,k)  + DeltaN;   % vypocet Xn+1

end     % norma posledniho vektoru nebude vypocitana

% norma posledniho F(Xn)
FXn = double(subs(F,{x y},Xn(:,N)')); % dosazeni Xn
rez(N+1) = norm(FXn, 2);

% Obrazky
fig1 = figure();
hold on
grid on
title({texlabel(f);texlabel(g)})
fimplicit(f)
fimplicit(g)
plot(Xn(1,:), Xn(2,:), 'k*')
hold off

fig2 = figure();
hold on
grid on
title('Rezidua')
plot(1:N+1, rez, 'ro')
hold off


%%%%%%%%%%%%%%%%%%%%%%
% % priklady ze seminare
% f=x^2-4*x+y^2;
% g=y-exp(-x)-1;
% A = [0; -2]; B = [2; 0]; C = [4; 1];
%
%
% f= x^2/9+y^2-1 ;
% g=y-log(x);
% A = [1.5; 1];B = [-1; 1]; C= [3; 1];
%
% f = 2*x - y*(1+x);
% g = x^2+5*y-1 ;
% A = [-1; 1];B = [0; 1];C = [0; 2];
%
% f= 1/ x-y^2;
% g= 2*x^2+y-4;
% A = [0; 1];
% B = [-0.5; -1];
% C = [1; 2];
%
% f= 1/(2*x)-y;
% g = x^2+4*y^2-4;
% X0=[1;0];