Nice! The Wikipedia article lists the solution that I call the "magical" one.
The trick to find it is to note that the AGM satisfies the functional equation
agm(x,y) = agm((x+y)/2,sqrt(xy))
To find a formula for the AGM, one would have to find a function that obeys such law. Once the answer has been given, it's simple. As the Wikipedia article shows, the elliptic integral is an invariant of the AGM iteration, so it should remain constant at every step. When the two numbers are the same, the elliptic integral is just the inverse. Therefore, once it's proven that the sequence converges, comparing the values of the elliptic integral at the beginning and "at infinity" one concludes that the AGM should be the inverse of the elliptic integral of the initial terms.
The brute force approach is like this. Observe that the AGM is homogeneous:
agm(kx,ky) = k*agm(x,y)
Because of this, the function you need to find actually depends on only one variable, not two. In principle, you can expand it in a Taylor series and use the functional equation to compare term by term. It turns out that you can prove this Taylor series is the one from the elliptic integral. The cleanest way I know of doing this expansion is by writing:
agm(1+x,1-x) = agm(1,(1-x)/(1+x))/(1+x) = agm(1,sqrt(1-4x/(1+x)²))/(1+x)
Now, agm(1+2sqrt(x)/(1+x),1-2sqrt(x)/(1+x)) = agm(1,sqrt(1-4x/(1+x)²)), so you can write the functional equation
(1+x)/agm(1+x,1-x) = 1/agm(1+2sqrt(x)/(1+x),1-2sqrt(x)/(1+x))
Now, you can expand 1/agm(1+x,1-x) in a Taylor series. It is simple because since it's an even function the odd exponents vanish. To expand the RHS you need the Taylor series for (1+x)
-n, which you can get by expanding 1/(1+x) as a geometric series and deriving the function. After some calculations you can get all coefficients and get the answer as a hypergeometric series. This series is exactly the elliptic integral and you could prove it using well known identities.
Elliptic functions are one of my favorite subjects because there are many different approaches to prove things. The variety of answers to this problem is one aspect of this. You can also recast the functional equation for the AGM as the hypergeometric differential equation, and then prove the answer by uniqueness of the IVP. There's even an approach using
theta functions. It turns out that the whole AGM sequence can be expressed as a ratio of squares of Jacobi theta functions, and you can infer all the properties from that!