Using the [definitions from Wikipedia](https://en.wikipedia.org/wiki/Lagrange_multiplier) 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`.