How to derive the hyperbola giving the foci and the fixed differene Given the two foci coordinates $(x_1,y_1)$ and $(x_2,y_2)$ of the hyperbola and the fixed difference distance, how can I derive its function to be able to draw it.

Let $f$ denote the distance between the foci, and let $d$ denote the constant difference between the distances from any point on the given hyperbola $H$ and the foci, which coincides with the distance between the vertices. The shape of the hyperbola (but not its size) is determined precisely by the ratio of these,
$$\epsilon := \frac{f}{d},$$
a parameter we call the eccentricity of the hyperbola. Now, one branch of an west-east facing hyperbola with eccentricity $\epsilon$ and distance $d$ between the vertices can be parameterized in polar coordinates (with the focus "inside" the parameterized branch at the origin) by
$$r(\theta) = \frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)},$$
which by construction is congruent $H$ via a rotation and a translation. Now, converting to rectangular coordinates gives the parameterization
$$\left\{\begin{array}{l} x(\theta) = \frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \cos \theta \\ y(\theta) = \frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \sin \theta\end{array}\right.$$
Next, we shift the hyperbola so that its center is at the origin. In this case, the distance to the focus is $\frac{1}{2}f$, so the shifted hyperbola is
$$\left\{\begin{array}{l} x(\theta) = \frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \cos \theta - \frac{f}{2}\\ y(\theta) = \frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \sin \theta\end{array}\right.$$
Now, to orient the hyperbola in the same direction, we need to apply a rotation transformation. Let $X := x_2 - x_1$ and $Y := y_2 - y_1$. We need to rotate the horizontal direction to the direction given by the vector $(X, Y)$ without scaling, and we can do that by applying the rotation
$$\frac{1}{X^2 + Y^2} \left(\begin{array}{cc} X & Y \\ -Y & X\end{array}\right).$$
Note that, by the Pythagorean theorem, $X^2 + Y^2 = f^2$. The result of this transformation is
$$\left\{\begin{array}{l} x(\theta) = \frac{1}{f^2}\left[X\left(\frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \cos \theta - \frac{f}{2}\right) + Y\left(\frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \sin \theta\right)\right]\\ y(\theta) = \frac{1}{f^2}\left[X\left(\frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \sin \theta\right) - Y\left(\frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \cos \theta - \frac{f}{2}\right)\right] \end{array}\right.$$
All that remains is to translate the hyperbola to coincide with $H$. We need to send the center $0$ of our hyperbola to the center of $H$, which we denote by $(\bar{x}, \bar{y})$, where $\bar{x} = \frac{1}{2}(x_1 + x_2)$ and $\bar{y} = \frac{1}{2}(y_1 + y_2)$, so the translation we want is $(x, y) \mapsto (x + \bar{x}, y + \bar{y})$. Composing with our parameterization gives
$$\left\{\begin{array}{l} x(\theta) = \frac{1}{f^2}\left[X\left(\frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \cos \theta - \frac{f}{2}\right) + Y\left(\frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \sin \theta\right)\right] + \bar{x}\\ y(\theta) = \frac{1}{f^2}\left[X\left(\frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \sin \theta\right) - Y\left(\frac{d(\epsilon^2 - 1)}{2(1 + \epsilon \cos \theta)} \cos \theta - \frac{f}{2}\right)\right] + \bar{y}\end{array}\right. .$$
It's possible to rearrange these steps some, and in particular I don't claim that this procedure is optimized for anything except (very arguably) geometric intuition.
