Let the points be $P_1(x_1, y_1)$ and $P_2(x_2, y_2)$, assumed to lie on an ellipse of semiaxes $a$ and $b$ with the $a$ axis making angle $\alpha$ to the $x$ axis.

Following @joriki, we rotate the points $P_i$ by $-\alpha$ into points
$$Q_i(x_i \cos(\alpha) + y_i \sin(\alpha), y_i \cos(\alpha) - x_i \sin(\alpha)).$$
We then rescale them by $(1/a, 1/b)$ to the points
$$R_i(\frac{x_i \cos(\alpha) + y_i \sin(\alpha)}{a}, \frac{y_i \cos(\alpha) - x_i \sin(\alpha)}{b}).$$

These operations convert the ellipse into a unit circle and the points form a chord of that circle. Let us now translate the midpoint of the chord to the origin: this is done by subtracting $(R_1 + R_2)/2$ (shown as $M$ in the figure) from each of $R_i$, giving points
$$S_1 = (R_1 - R_2)/2, \quad S_2 = (R_2 - R_1)/2 = -S_1$$
each of length $c$. Half the length of that chord is
$$c = ||(R_1 - R_2)||/2 = ||S_1|| = ||S_2||,$$
which by assumption lies between $0$ and $1$ inclusive. Set
$$s = \sqrt{1-c^2}.$$
The origin of the circle is found by rotating either of the $S_i$ by 90 degrees (in either direction) and rescaling by $s/c$, giving up to two valid solutions $O_1$ and $O_2$. (Rotation of a point $(u,v)$ by 90 degrees sends it either to $(-v,u)$ or $(v,-u)$.) For example, in the preceding figure it is evident that rotation $R_1$ by -90 degrees around $M$ and scaling it by $s/c$ will make it coincide with the circle's center. Reflecting the center about $M$ (which gives $2M$) produces the other possible solution.
Unwinding all this requires us to do the following to the $O_i$:
- Translate by $(R_1+R_2)/2$,
- Scale by $(a,b)$, and
- Rotate by $\alpha$.
The cases $c \gt 1$, $c = 1$, and $c=0$ have to be treated specially. The first gives no solution, the second a unique solution, and the third infinitely many.
FWIW, here's a Mathematica 7 function. The arguments p1 and p2 are length-2 lists of numbers (i.e., point coordinates) and the other arguments are numbers. It returns a list of the possible centers (or Null if there are infinitely many).
f[\[Alpha]_, a_, b_, p1_, p2_] := Module[
{
r, s, q1, q2, m, t, \[Gamma], u, r1, r2, x, v
},
(* Rotate to align the major axis with the x-axis. *)
r = RotationTransform[-\[Alpha]];
(* Rescale the ellipse to a unit circle. *)
s = ScalingTransform[{1/a, 1/b}];
{q1, q2} = s[r[#]] & /@ {p1, p2};
(* Compute the half-length of the chord. *)
\[Gamma] = Norm[q2 - q1]/2;
(* Take care of special cases. *)
If[\[Gamma] > 1, Return[{}]];
If[\[Gamma] == 0, Return[Null]];
If[\[Gamma] == 1,
Return[{InverseFunction[Composition[s, r]][(q1 + q2)/2]}]];
(* Place the origin between the two points. *)
t = TranslationTransform[-(q1 + q2)/2];
(* This ends the transformations.
The next steps find the centers. *)
(* Rotate the points 90 degrees. *)
u = RotationTransform [\[Pi]/2];
(* Rescale to obtain the possible centers. *)
v = ScalingTransform[{1, 1} Sqrt[1 - \[Gamma]^2]/\[Gamma]];
x = v[u[t[#]]] & /@ {q1, q2};
(* Back-transform the solutions. *)
InverseFunction[Composition[t, s, r]] /@ x
];