clear all, close all, clc;
p=@(x) 2*x.^2+1;
q=@(x) x.^4-x.^2;
z1=@(x) -p(x)/2+sqrt((p(x)/2).^2-q(x)); % Hade glömt (x) i sista q-termen
% z2=@(x) -p(x)/2-sqrt((p(x)/2).^2-q(x)); % Komplexa rötter
y11=@(x) sqrt(2*z1(x));
%y12=@(x) -sqrt(2*z1(x));
%y21=@(x) sqrt(2*z2(x));
%y22=@(x) -sqrt(2*z2(x));

x=-1:0.01:1;
y=y11(x);
plot(x,y,'b',x,-y,'b')
hold on;
y=x.^2-1/2;
plot(x,y)

% Kolla startgissningar ... (-0.3 , -0.3) & (-0.8 , 0.35)

F=@(X) [(X(1)^2+X(2)^2/2)^2-X(1)^2+X(2)^2/2;
    X(2)-X(1)^2+1/2];

J=@(X) [2*(X(1)^2+X(2)^2/2)*2*X(1)-2*X(1) 2*(X(1)^2+X(2)^2/2)*X(2)+X(2);
    -2*X(1) 1];

XX0=[-0.3 -0.8;-0.3 0.35];

for i=1:2
  X0=XX0(:,i);
  delta=inf;
  while norm(delta)> 1e-5
      f=F(X0);
      j=J(X0);
      delta=j\f;
      X0=X0-delta;
  end
  disp(['Rot ' num2str(i*2-1) ' x=' num2str(X0(1)) ' y=' num2str(X0(2))])
  disp(['Rot ' num2str(i*2) ' x=' num2str(-X0(1)) ' y=' num2str(X0(2))])
  plot(X0(1),X0(2),'*r',-X0(1),X0(2),'*r')
end