(* START of the Interior Point Method *)
x=0; y=1;s=1;z=1; μ=10; 
pts = {{x,y}};
A={{12(x-2)^2 + 2 - 2z,-4,0, -2x},{-4,8,0, 1},{0,0,-z, -s},{2x,-1,1, 0}};
b={-4(x-2)^3 - 2(x-2y)+2z*x,4(x-2y)-z,s*z+μ,-x^2+y-s};
step =LinearSolve[A,b];
x = x + step[[1]];y = y + step[[2]] ;s = s + step[[3]];z = z + step[[4]];
{N[x],N[y],N[s],N[z]} 
pts = Union[pts,{{x,y}}];

(* IPM iterations - run several times ... *)
μ=μ *0.1;
A={{12(x-2)^2 + 2 - 2z,-4,0, -2x},{-4,8,0, 1},{0,0,-z, -s},{2x,-1,1, 0}} ;
b={-4(x-2)^3 - 2(x-2y)+2z*x,4(x-2y)-z,s*z+μ,-x^2+y-s};
step =LinearSolve[A,b];
x = x + step[[1]];
y = y + step[[2]];
s = s + step[[3]];
z = z + step[[4]];
μ
{N[x],N[y],N[s],N[z]}
pts = Union[pts,{{x,y}}];

(* Use standard Mathematica function ... *)
Clear[x]; Clear[y];
FindMinimum[{(x-2)^4+(x-2y)^2,x^2-y ≤ 0},{{x,0},{y,1}}]

(* Plot IPM points *)
ContourPlot[(x-2)^4+(x-2y)^2,{x,0,1.5},{y,0,2},Epilog->{Red,PointSize[0.02],
{Thick,Line[pts]},Point[pts]},RegionFunction->Function[{x,y,z},x^2-y ≤ 0]]

(* Mathematica points *)
pts1 =Reap[FindMinimum[{(x-2)^4+(x-2y)^2,x^2-y ≤ 0},{{x,0},{y,1}},StepMonitor:>Sow[{x,y}]]][[2,1]];
pts1 =Join[{{0,1}},pts1];
ContourPlot[(x-2)^4+(x-2y)^2,{x,0,1.5},{y,0,2},Epilog->{Green,PointSize[0.02],{Thick,Line[pts1]},Point[pts1]},RegionFunction->Function[{x,y,z},x^2-y ≤ 0]]