There have been several responses to my posting on variation of parameters, so I decided to reply to the list as a whole.
First, and sadly, I have run out of reprints of the article I cited, so I'll give its gist below.
Lagrange did have a version of variation of parameters, but it bore little resemblance to what we know by that term today. That connection was made by Cauchy, and our use of the technique derives from Cauchy.
I have few details of the anticipation by Euler. It is mentioned by Enestrom in a footnote to Vessiot's article in the Encyclopedie (French, not German edition). See Enc. Sci. Math. 2(16), pp. 58ff, f/n 166. I have had no time since my last posting to add to this. The equation discussed was
and it arose in connection with a discussion of tides. This seems very close to modern usage, but the Euler source is an obscure one.
I was unaware of the references to Bernoulli and Leibniz. I'll try to follow them up -- when I get time! (A new term has just begun and we all have massive teaching loads!)
Lagrange introduced his technique in 1774 in a discussion of planetary orbits (Oeuvres 4, 109-147) but expanded on it in 1781 (Oeuvres 5, 123ff). For Cauchy, see Oeuvres Completes 2 ser. 6, 252-255.
Here, from a later text (Vol 3 of Price's 1856 Treatise on Infinitesimal Calculus), is an example of Lagrange's approach, which is essentially a perturbation technique.
Consider the equation of resisted fall
and the related unperturbed equation
The latter has the solution
so for the former set
In the case of the simpler equation we have
and so for the perturbed equation, we expect
For consistency impose the condition
Putting this all together
and I leave it to you all to complete the work and reach the full solution.
As I said, the roots of the method are present in this example (as in Lagrange's more complicated ones) but the routine technique we use today was extracted from it by Cauchy.
Which is not to say that Euler, Bernoulli and Leibniz didn't also independently and earlier entertain similar thoughts!