Using the definitions from Wikipedia rather literaly:

Clear[f, g]

f[x_, y_] := 6 x + 8 y
g[x_, y_] := x^2 + y^2 - 25

Solve[
 D[f[x, y] - lambda (g[x, y] - 0), {{x, y, lambda}}] == 0,
 {x, y}, lambda
]

(* Out: {{x -> -3, y -> -4}, {x -> 3, y -> 4}} *)

Using the definitions from Wikipedia rather literally:

Clear[f, g]

f[x_, y_] := 6 x + 8 y
g[x_, y_] := x^2 + y^2 - 25

Solve[
 D[f[x, y] - lambda g[x, y], {{x, y, lambda}}] == 0,
 {x, y}, lambda
]

(* Out: {{x -> -3, y -> -4}, {x -> 3, y -> 4}} *)

Unpacking the code, we start with defining the function to be optimized, and the constraint:

Clear[f, g]

f[x_, y_] := 6 x + 8 y
g[x_, y_] := x^2 + y^2 - 25

The method of Lagrange multipliers assures us that the extrema of the original function $f$ are stationary points for $\Lambda = f-\lambda \ g$. In order to find those, we calculate the gradient of this auxiliary function $\Lambda$, and set it to zero:

D[f[x, y] - lambda g[x, y], {{x, y, lambda}}] == 0

For a function $f$ of $n$ variables, this gives us a set of $n+1$ equations in $n+1$ variables, the original ones plus the multiplier, which we can solve using Solve.