If I have a program that can compute FFT's for sizes that are powers of 2, how can I use it to compute FFT's for other sizes?

I have read that I can supposedly zero-pad the original points, but I'm lost as to how this gives me the correct answer. If I transform a function with 5 points, I expect to get back a 5-point DFT, and I don't understand how I can extract that from the 8-point DFT (which was calculated with zero-padding).

Sure, you can use a radix-2 FFT to compute FFTs for lengths not a power of 2 (but it is not as efficient as using methods specifically tailored to the factors of the sequence length). However, it's not as easy as merely padding the original array. The right way of going about it, due to Rabiner, Schafer, and Rader, is termed the "chirp-z transform". (This is a slightly condensed version of the paper previously linked to.)
To motivate the chirp-z idea, consider the DFT
$$X_k=\sum_{n=0}^{N-1} x_n \left(e^{-2\pi i/N}\right)^{k n},\qquad k = 0,\dots,N-1$$
The key step, attributed to Bluestein, considers the algebraic identity
$$kn=\frac12(k^2+n^2-(k-n)^2)$$
Substituting into the DFT and expanding out the powers yields
$$\begin{align*}&\sum_{n=0}^{N-1} x_n \left(e^{-2\pi i/N}\right)^{\frac12(k^2+n^2-(k-n)^2)}\\=&\left(\left(e^{-2\pi i/N}\right)^{k^2/2}\right)\sum_{n=0}^{N-1} \left(x_n \left(e^{-2\pi i/N}\right)^{n^2/2}\left(e^{-2\pi i/N}\right)^{-\frac{(k-n)^2}2}\right)\end{align*}$$
and we then recognize that the summation expression is precisely in the form of a convolution. The factor $\left(e^{-2\pi i/N}\right)^{k^2/2}$ is what is termed as a "chirp"; hence the name "chirp-z transform".
Chirp-z thus consists of three stages:

taking the Hadamard (componentwise) product of the original sequence with the chirp
convolving with the reciprocal chirp (which of course is equivalent to the inverse FFT of the product of the FFTs of the Hadamard product and the reciprocal chirp)
taking the Hadamard product of that result with a chirp. We see that three or so FFTs are needed (barring the possibility of an exploitable symmetry being present).

MATLAB's Signal Processing toolbox implements this algorithm as the function czt(). See the code in the corresponding M-file for implementation details.

The Code
I (the original poster) ended up making a Python version using SciPy, so I thought I'd edit this answer to include the code:

from scipy import *

def czt(x, m=None, w=None, a=None):
    # Translated from GNU Octave's czt.m

    n = len(x)
    if m is None: m = n
    if w is None: w = exp(-2j * pi / m)
    if a is None: a = 1

    chirp = w ** (arange(1 - n, max(m, n)) ** 2 / 2.0)
    N2 = int(2 ** ceil(log2(m + n - 1)))  # next power of 2
    xp = append(x * a ** -arange(n) * chirp[n - 1 : n + n - 1], zeros(N2 - n))
    ichirpp = append(1 / chirp[: m + n - 1], zeros(N2 - (m + n - 1)))
    r = ifft(fft(xp) * fft(ichirpp))
    return r[n - 1 : m + n - 1] * chirp[n - 1 : m + n - 1]

@vectorize  # Rounds complex numbers
def cround(z, d=None): return round(z.real, d) + 1j * round(z.imag, d)

arr = [1.0, 2.0, 3.0]

print cround(czt(arr), 4)       # [ 6.0+0.j    -1.5+0.866j -1.5-0.866j]
print cround(fft(arr), 4)       # [ 6.0+0.j    -1.5+0.866j -1.5-0.866j]

