Finding a possible Lagrangian corresponding to a differential equation

A typical problem in the calculus of variations is to extremize a functional

$$ J[y]=int f(x,y,y’) , mathrm{d} x.$$

This usually involves solving the Euler-Lagrange equation

$$frac{mathrm{d}}{mathrm{d}x} frac{partial f}{partial y’}-frac{partial f}{partial y}=0.$$

The inverse of this problem is given a differential equation, determine $f$ such that a function $y$ is a solution to the original differential equation if and only if $y$ is a solution to the Euler-Lagrange equation. For some cases, like second-order linear differential equations, there are known solutions to this problem. In fact, Wolfram|Alpha can solve this exact problem. For example, running

WolframAlpha["y''+2y'-y=0", {{"PossibleLagrangian", 1}, "ComputableData"}]

gives one such solution:

Hold[ℒ[y', y, x] == 1/2 (E^(2 x) y^2 + E^(2 x) (y')^2)]

Is there a way to do this directly in Mathematica, without using Wolfram|Alpha or resorting to known formulas?

I’ve written a short function based on a known formula (The Calculus of Variations by Brunt, section 3.4), which gives the same answer as WA:

(* y'' + P y' + Q y - G == 0 *)

PossibleLagrangian[P_, Q_, G_, x_] := 
 Block[{p = Exp[Integrate[P[t], {t, 0, x}]], q = Q[x] p, g = G[x] p},
  1/2 (p (y')^2 - q y^2 + 2 g y)]

(* y'' + 2y' - y == 0 *)
PossibleLagrangian[2 &, -1 &, 0 &, x]
1/2 (E^(2 x) y^2 + E^(2 x) (y')^2)

 Names["*`*Lagrang*"]
(*{"NDSolve`FEM`HangingNodeLagrangeMultipliers", 
   "NDSolve`FEM`LagrangeMultipliers", "Region`LagrangeNearest"}*)