function fun
format long
intvalinit('DisplayInfsup')
%General Case -- DZ1
disp('General Case -- DZ1')
syms a b c d
f1=a^4-b*c*d;
f2=b^4-a*c*d;
f3=c^4-a*b*d;
f4=d^4-a*b*c;
var=[a,b,c,d];
pstar=rand(1,4)*10^(-4);
sys=[f1,f2,f3,f4];

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%General Case -- DZ2
disp('General Case -- DZ2')
syms x y z
f1=x^4;
f2=x^2*y+y^4;
f3=z+z^2-7*x^3-8*x^2;
sys=[f1,f2,f3];
var=[x,y,z];
pstar=rand(1,3)*10^(-4);
pstar(3) = pstar(3)-1;

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- DZ3
disp('General Case -- DZ3')
syms x y
f1=14*x+33*y-3*sqrt(5)*x^2-12*sqrt(5)*x*y-12*sqrt(5)*y^2-6*sqrt(5)+x^3+6*x^2*y+12*x*y^2+8*y^3+sqrt(7);
f2=41*x-18*y-sqrt(5)+8*x^3-12*x^2*y+6*x*y^2-y^3+12*sqrt(7)*x*y-12*sqrt(7)*x^2-3*sqrt(7)*y^2-6*sqrt(7);
var = [x,y];
pstar = [2*sqrt(7)/5+sqrt(5)/5,-sqrt(7)/5+2*sqrt(5)/5] + rand(1,2)*10^(-4);
sys =[f1,f2];

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))


%General Case -- cbms1
disp('General Case -- cbms1')
syms x y z
f1=x^3-y*z;
f2=y^3-x*z;
f3=z^3-x*y;
sys=[f1,f2,f3];
var=[x,y,z];
pstar=rand(1,3)*10^(-4);

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%General Case -- cbms2
disp('General Case -- cbms2')
syms x y z
f1=x^3-3*x^2*y+3*x*y^2-y^3-z^2;
f2=y^3-3*y^2*z+3*y*z^2-z^3-x^2;
f3=z^3-3*z^2*x+3*z*x^2-x^3-y^2;
sys=[f1,f2,f3];
var=[x,y,z];
pstar=rand(1,3)*10^(-4);

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%General Case -- mth191
disp('General Case -- mth191')
syms x y z
f1=x^3+y^2+z^2-1;
f2=x^2+y^3+z^2-1;
f3=x^2+y^2+z^3-1;
sys=[f1,f2,f3];
var=[x,y,z];
pstar=rand(1,3)*10^(-4);
pstar(2) = pstar(2)+1;

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%General Case -- KSS
disp('General Case -- KSS')
syms x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
f1=x1^2-x1+x2+x3+x4+x5-4;
f2=x2^2+x1-x2+x3+x4+x5-4;
f3=x3^2+x1+x2-x3+x4+x5-4;
f4=x4^2+x1+x2+x3-x4+x5-4;
f5=x5^2+x1+x2+x3+x4-x5-4;
f6=x6^2+x1+x2+x3+x4+x5-2*x6-4;
f7=x7^2+x1+x2+x3+x4+x5-2*x7-4;
f8=x8^2+x1+x2+x3+x4+x5-2*x8-4;
f9=x9^2+x1+x2+x3+x4+x5-2*x9-4;
f10=x10^2+x1+x2+x3+x4+x5-2*x10-4;


sys=[f1,f2,f3,f4,f5,f6,f7,f8,f9,f10];
var=[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10];
pstar=rand(1,10)*10^(-4)+1;


tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- RuGr09
disp('Special Case -- RuGr09')
syms a b
f1=a^2*b-a*b^2;
f2=a-b^2;

var = [a,b];
pstar = rand(1,2)*10^(-4);

sys =[f1,f2];

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- LZ
disp('Special Case -- LZ')
sys=[myseq('x[i]^2+x[i]','i',1:99)+myseq('-x[i]','i',2:100);myseq('x[j]^3','j',100)].';
var=[myseq('x[i]','i',1:100)].';
pstar = rand(1,100)*10^(-4);

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- Ojika1
disp('Special Case -- Ojika1')
syms x y
f1=x^2+y-3;
f2=x+1/8*y^2-3/2;

var = [x,y];
pstar = rand(1,2)*10^(-4) + [1,2];

sys =[f1,f2];

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- Ojika2
disp('Special Case -- Ojika2')
syms x y z
f1=x^2+y+z-1;
f2=x+y^2+z-1;
f3=x+y+z^2-1;

var = [x,y,z];
pstar= rand(1,3)*10^(-4) + [1,0,0];

sys =[f1,f2,f3];
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- Ojika3
disp('Special Case -- Ojika3')
syms x y z
f1=x+y+z-1;
f2=2*x^3+5*y^2-10*z+5*z^2+5;
f3=2*x+2*y+z^2-1;

var = [x,y,z];
pstar = rand(1,3)*10^(-4) + [-2.5,2.5,1];

sys =[f1,f2,f3];
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- Ojika4
disp('Special Case -- Ojika4')
syms x y z
f1=x+z*x^3+x*z*y^2-x*z;
f2=10*y-2*y*z*x^2-z*y^3-y*z;
f3=-6*z^2*x^4-3*x^2*y^2*z^2-z^2*x^2+28*z*x^2-3*z^2*y^4+2*z^2*y^2+7*z*y^2+z^2-11*z+10;

var = [x,y,z];
pstar = rand(1,3)*10^(-4) + [0,0,1];

sys =[f1,f2,f3];

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- Decker2
disp('Special Case -- Decker2')
syms x y
f1=x+y^3;
f2=x^2*y-y^4;
var = [x,y];
pstar = rand(1,2)*10^(-4);

sys =[f1,f2];

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%General Case -- Caprasse
disp('General Case -- Caprasse')
syms x1 x2 x3 x4
f1=-x1^3*x3+4*x1*x2^2*x3+4*x1^2*x2*x4+2*x2^3*x4+4*x1^2-10*x2^2+4*x1*x3-10*x2*x4+2;
f2=-x1*x3^3+4*x2*x3^2*x4+4*x1*x3*x4^2+2*x2*x4^3+4*x1*x3+4*x3^2-10*x2*x4-10*x4^2+2;
f3=x2^2*x3+2*x1*x2*x4-2*x1-x3;
f4=2*x2*x3*x4+x1*x4^2-x1-2*x3;
sys=[f1,f2,f3,f4];
var=[x1,x2,x3,x4];
pstar=rand(1,4)*10^(-4)+rand(1,4)*10^(-4)*sqrt(-1)+[2,-sqrt(-3),2,sqrt(-3)];

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- cyclic9 
disp('Special Case -- cyclic9 ')
syms y1 y2 y3 y4 y5 y6 y7 y8 y9
f1=y1+y2+y3+y4+y5+y6+y7+y8+y9;
f2=y1*y2+y2*y3+y3*y4+y4*y5+y5*y6+y6*y7+y7*y8+y8*y9+y9*y1;
f3=y1*y2*y3+y2*y3*y4+y3*y4*y5+y4*y5*y6+y5*y6*y7+y6*y7*y8+y7*y8*y9+y8*y9*y1+y9*y1*y2;
f4=y1*y2*y3*y4+y2*y3*y4*y5+y3*y4*y5*y6+y4*y5*y6*y7+y5*y6*y7*y8+y6*y7*y8*y9+y7*y8*y9*y1+y8*y9*y1*y2+y9*y1*y2*y3;
f5=y1*y2*y3*y4*y5+y2*y3*y4*y5*y6+y3*y4*y5*y6*y7+y4*y5*y6*y7*y8+y5*y6*y7*y8*y9+y6*y7*y8*y9*y1+y7*y8*y9*y1*y2+y8*y9*y1*y2*y3+y9*y1*y2*y3*y4;
f6=y1*y2*y3*y4*y5*y6+y2*y3*y4*y5*y6*y7+y3*y4*y5*y6*y7*y8+y4*y5*y6*y7*y8*y9+y5*y6*y7*y8*y9*y1+y6*y7*y8*y9*y1*y2+y7*y8*y9*y1*y2*y3+y8*y9*y1*y2*y3*y4+y9*y1*y2*y3*y4*y5;
f7=y1*y2*y3*y4*y5*y6*y7+y2*y3*y4*y5*y6*y7*y8+y3*y4*y5*y6*y7*y8*y9+y4*y5*y6*y7*y8*y9*y1+y5*y6*y7*y8*y9*y1*y2+y6*y7*y8*y9*y1*y2*y3+y7*y8*y9*y1*y2*y3*y4+y8*y9*y1*y2*y3*y4*y5+y9*y1*y2*y3*y4*y5*y6;
f8=y1*y2*y3*y4*y5*y6*y7*y8+y2*y3*y4*y5*y6*y7*y8*y9+y3*y4*y5*y6*y7*y8*y9*y1+y4*y5*y6*y7*y8*y9*y1*y2+y5*y6*y7*y8*y9*y1*y2*y3+y6*y7*y8*y9*y1*y2*y3*y4+y7*y8*y9*y1*y2*y3*y4*y5+y8*y9*y1*y2*y3*y4*y5*y6+y9*y1*y2*y3*y4*y5*y6*y7;
f9=y1*y2*y3*y4*y5*y6*y7*y8*y9-1;
sys=[f1,f2,f3,f4,f5,f6,f7,f8,f9];
var=[y1,y2,y3,y4,y5,y6,y7,y8,y9];
pstar=[-.9396926-.3420201*sqrt(-1),-2.4601472-.8954204*sqrt(-1), -.3589306-.1306401*sqrt(-1),-.9396926-.3420201*sqrt(-1),.3589306+.1306401*sqrt(-1),2.4601472+.8954204*sqrt(-1),-.9396926-.3420201*sqrt(-1),.3589306+.1306401*sqrt(-1), 2.4601472+.8954204*sqrt(-1)];
appsol=zeros.';
tol= 0.005;

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Special Case -- LZ500
disp('Special Case -- LZ500')
sys=[myseq('x[i]^2+x[i]','i',1:499)+myseq('-x[i]','i',2:500);myseq('x[j]^3','j',500)].';
var=[myseq('x[i]','i',1:500)].';
pstar = rand(1,500)*10^(-4);

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- LZ1000
disp('Special Case -- LZ1000')
sys=[myseq('x[i]^2+x[i]','i',1:999)+myseq('-x[i]','i',2:1000);myseq('x[j]^3','j',1000)].';
var=[myseq('x[i]','i',1:1000)].';
pstar = rand(1,1000)*10^(-4);

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%Special Case -- LZ2000
disp('Special Case -- LZ2000')
sys=[myseq('x[i]^2+x[i]','i',1:1999)+myseq('-x[i]','i',2:2000);myseq('x[j]^3','j',2000)].';
var=[myseq('x[i]','i',1:2000)].';
pstar = rand(1,2000)*10^(-4);

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,10^(-6),0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))
%simple1
disp('simple1')
syms x y z w t
f1=y+x^3*y+z^5;
f2=-x+x^2*y;
f3=3*x^2*y*w+3*x^2*y+2*x*y*t-t;
f4=x^3*w+x^3+x^2*t+2*y+w+1;
f5=z^4*w+z^4;
sys=[f1,f2,f3,f4,f5];
var=[x,y,z,w,t];
pstar=rand(1,5)*10^(-4)+[0,0,0,-1,0];
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%simple2
disp('simple2')
 syms a b c d
f1=a+c^2;
f2=b;
f3=c+d^2;
f4=d^5;
sys=[f1,f2,f3,f4];
var=[a,b,c,d];
pstar=rand(1,4)*10^(-4);
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%simple3
disp('simple3')
syms a b
f1=a+b;
f2=a^10+b^10;
sys=[f1,f2];
var=[a,b];
pstar=rand(1,2)*10^(-4);
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%simple4
disp('simple4')
syms x1 x2 x3 x4 x5 x6
f1=3*x1+x2*x4-x6^4;
f2=-3*x1*x4-3*x3*x5+x3*x6-x4*x5-3*x2-2*x6;
f3=-x1*x3-2*x6^2-3*x3+x4*x5*x6^4;
f4=-3*x1^3-x1*x2-2*x2*x5^3-x6^2+2*x4+x1;
f5=-3*x1*x4^2-3*x1*x6-2*x5*x6-3*x5;
f6=-2*x1*x5^2-x2^2*x4+x3^3+3*x3^2*x6;
sys=[f1,f2,f3,f4,f5,f6];
var=[x1,x2,x3,x4,x5,x6];
pstar=rand(1,6)*10^(-4);
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('cluster1')
syms x y 
f1=10^6*x^2;
f2=10^6*y^2;
sys=[f1,f2];
var=[x,y];
pstar = rand(1,2)*10^(-4);

tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,1);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%%%%%%%%%%%%%%%%%%%%%%
disp('cluster2')
syms x y 
f1=10^6*x^2;
f2=y^2;
sys=[f1,f2];
var=[x,y];
pstar = rand(1,2)*10^(-4);
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,1);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%%%%%%%%%%%%%%
disp('cluster3')
syms x y 
f1=10^6*x^2+x;
f2=10^6*x^2+y;
sys=[f1,f2];
var=[x,y];
pstar = rand(1,2)*10^(-4);
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,0);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))


%%%%%%%%%%%%%%%%%

disp('cluster4')
syms x y 
f1=10^6*x^2+y;
f2=y;
sys=[f1,f2];
var=[x,y];
pstar = rand(1,2)*10^(-4);
tic;[X,xs,g,G,newvar]=refine(sys,var,pstar,0.01,1);toc
m=max(size(xs))
n=max(size(var));
p=xs';
if ~isempty(G)
Gstar=subs(G,newvar(n+1:m),p(n+1:m));
sys=[sys,Gstar];
end
err=maxerr(sys,var,p(1:n))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

