Catenary Between Two Points (Hanging Chain) This is one of several challenges left for the community by Calvin's Hobbies.

The curve that an idealised hanging rope or chain makes is a catenary.  


Image by Bin im Garten, via Wikimedia Commons. Used under the CC-By-SA 3.0 licence.

Write a program that will draw a catenary, as an image, in quadrant 1 of the plane given two points (x1,y1), (x2,y2), and the "rope length" L. L will be greater than the distance between the two points.

You must also draw axes on the left and bottom sides of the image (400x400 px min) for scale. Only draw the quadrant from x and y in range 0 to 100. (You may assume the points are in range.)

Dots, or something similiar, should be drawn at the (x1,y1), (x2,y2) endpoints to distinguish them. The curve should only be drawn in the space between these points.
Python + NumPy + Matplotlib, 1131

Just to get us started, here's an attempt that uses no knowledge of calculus or physics other than the fact that the catenary minimizes the energy of a chain. Hey, my algorithm may not be efficient, but at least it's not implemented efficiently either!

import math
import random
import numpy as np
import matplotlib.pyplot as plt
length, x0, y0, x1, y1 = input(), input(), input(), input(), input()
chain = np.array([[x0] + [length / 1000.]*1000, [y0] + [0.] * 1000])
def rotate(angle, x, y):
 return x * math.cos(angle) + y * math.sin(angle), -x * math.sin(angle) + y  * math.cos(angle)
def eval(chain, x1, y1):
 mysum = chain.cumsum(1)
 springpotential = 1000 * ((mysum[0][-1] - x1) ** 2 + (mysum[1][-1] - y1)  ** 2)
 potential = mysum.cumsum(1)[1][-1]
 return springpotential + potential
def jiggle(chain, x1, y1):
 for _ in xrange(100000):
  pre = eval(chain, x1, y1)
  angle = random.random() * 2 * math.pi
  index = random.randint(1,1000)
  chain[0][index], chain[1][index] = rotate(angle, chain[0][index], chain[1][index])
  if( pre &lt; eval(chain, x1, y1)):
   chain[0][index], chain[1][index] = rotate(-angle, chain[0][index], chain[1][index])
jiggle(chain, x1, y1)
sum = chain.cumsum(1)
x1 = 2 * x1 - sum[0][-1]
y1 = 2 * y1 - sum[1][-1]
jiggle(chain, x1, y1)
sum = chain.cumsum(1)
plt.plot(sum[0][1:], sum[1][1:])
plt.show()