(* We find the maximum of f(x,y,z) =  2x+y+z on the intersection of the plain x-y+2z = 1 and the paraboloid z = x^2+2y^2 *)
f=2x+y+z
g=x-y+2 z
h=x^2+2 y^2-z
ContourPlot3D[f,{x,-1,1},{y,-1,1},{z,-1,1}]
ContourPlot3D[{g==1,h==0,},{x,-2,1.5},{y,-1,1},{z,0,2}]

Gf=Grad[f,{x,y,z}]
Gg=Grad[g,{x,y,z}]
Gh=Grad[h,{x,y,z}]

sol=Solve[{Gf==alpha Gg+beta Gh,g==1,h==0}]

s1={x,y,z,f}/.sol[[1]] 
s1 //N

s2={x,y,z,f}/.sol[[2]] 
s2 //N

(* Approximate global extrema* )
(* Global maximum of 1.728 at (0.379,0.439,0.530 *)
(* Global minimum of -1.103 at (-0.879,-0.189,0.844) *)