function p = inverseinter2Droot(f,x1,x2,x3)
%Inverse quadratic Lagrange erpolation
p=1;
maxit=100,iter=0;
tol=1.0e-10;
es=0.00000001;
ea=0.001;
f1=f(x1);
f2=f(x2);
f3=f(x3);
if f1==0 p=x1;
elseif f2==0 p=x2;
elseif f3==0 p=x3;end 
if f1*f3>0 fprintf('Root is not existed in the given region');end
while (ea>es) && (iter<maxit) && f1~=0 && f2~=0 && f3~=0 
    p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
    fp=f(p);
    ea=abs(x3-p);
	if abs(f3)<tol  p=x3; end       
    if p>x3  x1=x3;f1=f3;x3=p;f3=fp;
    elseif(p<x3) x2=x3;f2=f3;x3=p;f3=fp;end
    iter=iter+1;
    end 
if(iter>=maxit) fprintf('Warning : MAximum number of iteration is exceeded \nresult may not be valid');end
end
