
function er=testNmultzero(h=0.1,p=2)
%we test local newton convergence for 
%f=@(x) ((x-2)^p)*(exp(x)+1) (default p=2)
%starting from x0=2+h (2 is a root of multiplicity p)
%INPUT
%h - starting error i.e. x0-2=h
%p - multiplicity of the root 2
%OUTPUT
%er - vector with errors .e .er(k) = x_k-2
%EXAMPLE
%double zero
%er=testNmultzero(h=0.1,p=2); n=length(er);er(2:n)./er(1:n-1)
%single zero
%er=testNmultzero(h=0.1,p=1); n=length(er);er(2:n)./er(1:n-1)
%triple zero
%er=testNmultzero(h=0.1,p=3); n=length(er);er(2:n)./er(1:n-1)
%noninteger multiplicity p=2.5
%er=testNmultzero(h=0.1,p=2.5); n=length(er);[er(2:n)./er(1:n-1);1-1/p]
f=@(x) ((x-2)^p)*(exp(x)+1) ;
df=@(x) (p*(x-2)^(p-1)).*(exp(x)+1) + ((x-2)^p).*exp(x);

x=2+h;
maxi=20;it=0;
fx=f(x);
err=h;
er=zeros(maxi,1);
printf("[it] err/ep err/(ep^2)  err/(ep^3) \n");
while((abs(fx)>1e-10) && (it<maxi))
x=x-fx/df(x);
ep=err;
er(++it)=err=x-2;
printf("[%d] %g %g %g \n",it,err/ep,err/(ep*ep),err/(ep*ep*ep));
fx=f(x);
end
er=er(1:it);
end