Python + NumPy + Matplotlib, 1131 203

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 < 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()

Edit: Since this is getting comments, lets do a 10 year anniversary re-golf!

from matplotlib.pyplot import*
def j(l,z,q,*x):m=999;c=np.array([[z,q]]*m);exec(("i=1+np.random.randint(m-1);v=x-c.sum(0)-[0,(m-i)/m];c[i]=v/m*l/v.dot(v)**.5;"*m*19+"x+=v;")*2);plot(*c.cumsum(0).T);show()

example usage: j(10, 1, 0., 4., 1.)

Instead of optimizing by picking random angles and hillclimbing the energy, this uses the closed form solution for the greedy update on a single link.

BBC Basic, 300 ASCII characters, tokenised filesize 260

  INPUTr,s,u,v,l:r*=8s*=8u*=8v*=8l*=8z=0REPEATz+=1E-3UNTILFNs(z)/z>=SQR(l^2-(v-s)^2)/(u-r)a=(u-r)/2/z
  p=(r+u-a*LN((l+v-s)/(l-v+s)))/2q=(v+s-l*FNc(z)/FNs(z))/2MOVE800,0DRAW0,0DRAW0,800CIRCLEu,v,8CIRCLEr,s,8FORx=r TOu
    DRAW x,a*FNc((x-p)/a)+q
  NEXT
  DEFFNs(t)=(EXP(t)-EXP(-t))/2
  DEFFNc(t)=(EXP(t)+EXP(-t))/2

Emulator at http://www.bbcbasic.co.uk/bbcwin/bbcwin.html

This has obviously been solved before, so the first thing I did was look what others have done.

The equation of a catenary centred at the origin is simply y=a*cosh(x/a). It becomes slightly more complicated if it is not centred at the origin.

Various sources say that if the length and endpoints are known the value for a must be determined numerically. There is an unspecified parameter h in the wikipedia article. So I found another site and basically followed the method here: http://www.math.niu.edu/~rusin/known-math/99_incoming/catenary

BBC Basic does not have sinh and cosh built in, so I defined two functions at the end of the program to calculate them using EXP

coordinates for lefthand point must be supplied before righthand point, OP confirmed this is OK. Length is given last. Values can be separated by commas or newlines.

Ungolfed code

  INPUT r,s,u,v,l

  REM convert input in range 0-100 to graphic coordinates in range 0-800 
  r*=8 s*=8 u*=8 v*=8 l*=8

  REM solve for z numerically
  z=0
  REPEAT
    z+=1E-3
  UNTIL FNs(z)/z>=SQR(l^2-(v-s)^2)/(u-r)

  REM calculate the curve parameters
  a=(u-r)/2/z
  p=(r+u-a*LN((l+v-s)/(l-v+s)))/2
  q=(v+s-l*FNc(z)/FNs(z))/2

  REM draw axes, 800 graphics units long = 400 pixels long (2 graphics units per pixel)
  MOVE 800,0
  DRAW 0,0
  DRAW 0,800

  REM draw markers at end and beginning of curve (beginning last, so that cursor is in right place for next step)
  CIRCLE u,v,8
  CIRCLE r,s,8

  REM draw curve from beginning to end
  FORx=r TOu
    DRAW x,a*FNc((x-p)/a)+q
  NEXT

  REM definitions of sinh and cosh
  DEF FNs(t)=(EXP(t)-EXP(-t))/2
  DEF FNc(t)=(EXP(t)+EXP(-t))/2

Python 2.7 + matplotlib, 424

Run as

python thisscript.py [length] [x0] [y0] [x1] [y1]

If I can assume that x0 is always smaller than x1 character count reduces to 398

from numpy import *
from pylab import *
from scipy.optimize import *
import sys
c=cosh
l,p,q,u,w=map(float,sys.argv[1:])
if p>u:
 p,q,u,w=u,w,p,q
h=u-p
v=w-q
a=brentq(lambda a:(2.*h/a*sinh(0.5*a))**2-l**2-v**2,1e-20,600)
b=brentq(lambda b:c(a*(1.-b))-c(a*b)-a*v/h,-600/a,600/a)
r=linspace(p,u,100)
plot([p,u],[q,w],'ro')
plot(r,h/a*c(((r-p)/h-b)*a)-h/a*c(a*b)+q,'k-')
gca().set_xlim((0,100))
gca().set_ylim((0,100))
show()

The magic number 600 you see appearing in some places is due to the fact that cosh(x) and sinh(x) start overflowing around x=710 (so 600 to keep some margin)

Basically I solve the problem in the frame where the catenary goes trough (0,0) and (x1-x0,(y1-y0)/(x1-x0)) and then remap to the original frame. This improves the numerical stability a lot.