Simpson's Rule like Method for Arc Length Approximation Simpson's Rule is of course a very important method for Numerical Integration.  The basic idea behind it is to the intervals as a quadratic curve, and to calculate the area under that curve (yes, I realize I'm simplifying it:  I don't think I need to go over the details).
What I'm curious is, if I want to Numerically calculate the arc length of a function over a given interval, is there a similar trick we can use to calculate it?  I'm aware of the basic method behind calculating arclength 
$\sum_{i=0}^n \sqrt(dt^2 + (f(i + 1)-f(i))^2) $
, but can't find anything helpful similar to the Simpson's rule.
Also, the only function we can use are those from the original function.
So for instance, if I want to find
$Arc Length = \sqrt(x) $ from x = 0 to 100, the only functions we can use values from is $f(x) = \sqrt(x) $
I bring this up because I've seen several sites that just say to use Simpson's Rule Integration with
$ \int_0^{100} \sqrt(1+ [f'(x)]^2)dx $
but that is NOT what I am referring to (that's cheating :)
Thanks in advance.
 A: Simpson's Rule is the result of Richardson's extrapolation applied to trapezoidal rule, as seen here. So, we can obtain its natural analog for the arclength by applying the same extrapolation to the sum-of-chordarcs formula. It goes like this: 
Let $n$ be an even integer. Take the sum of chords over $n$ chords:
$$L \approx \sum_{i=1}^{n} \sqrt{(x_{i }-x_{i-1})^2 + (f(x_{i })-f(x_{i-1}))^2 } \tag{1}$$
where $x_i = a+i (b-a)/n$, $i=0,\dots,n$. Then do the same for $n/2$ chords:
$$L \approx \sum_{i=1}^{n/2} \sqrt{(x_{2i}-x_{2i-2})^2 + (f(x_{2i})-f(x_{2i-2}))^2 } \tag{2}$$
Since the sum-of-chords formula has error of order $1/n^2$ (error comes from second derivative of $f$), the error in (2) should be about $4$ times the error in (1). So, divide (2) by $4$ and subtract it from (1). Result: 
$$  L \approx \frac43 \sum_{i=1}^{n} \sqrt{(x_{i }-x_{i-1})^2 + (f(x_{i })-f(x_{i-1}))^2 } \\ - 
\frac13 \sum_{i=1}^{n/2} \sqrt{(x_{2i}-x_{2i-2})^2 + (f(x_{2i})-f(x_{2i-2}))^2 }  \tag{3}$$
This improves precision in the same manner as Simpson's rule improves over trapezoidal method.

Example: The length of parabola $y=x^2$ between $(0,0)$ and $(10,100)$; the same graph as in your post, turned sideways. Using $n=100$ in (1),
I get within $8.3\cdot 10^{-4}$ of the exact length. Using the same $n=100$ in (3), I get 
within $10^{-9}$ of the exact length.
Huge improvement, with no extra evaluations of $f$ (the second sum in (3) uses the values of $f$ already found for the first sum). It appears that (3) is of 4th order of accuracy when $f$ is smooth, just as Simpson's rule.
The improvement is far less spectacular for $\sqrt{x}$ on $[0,100]$, which is not smooth at $0$. This interval calls for an adaptive method, because on most of it the function is almost flat but near $0$ its derivative blows up. This is true for both area and length computations. To make (3) adaptive, subdivide $[a,b]$ recursively, stopping 
when the results of (1) and (2) are sufficiently close to each other. Apply (3) on each subinterval you obtained, and sum the result.

The algorithm proposed by   Vincent  and Forsey, pointed out by Biswajit Banerjee in a comment, is likely to outperform (3) - after all,  they spent time thinking about particulars of length computation, while (3) is a straightforward extrapolation of sum-of-chords formula.
