If $|x_1 - x_2| > |y_1 - y_2|$, then the line will intersect the square in one of the sides, left or right. Then if $x_1 > x_2$, the point of intersection is on the left, so the point will be some $(x,y)$ with $x = x_2 - 15$. Since it is on the line $(x_2,y_2) + \lambda (x_2 - x_1, y_2 - y_1)$ we know that for some $\lambda$, $x_2 - 15 = x_2 + \lambda(x_2 - x_1)$ and $y = y_2 + \lambda(y_2 - y_1)$. The first equality gives you $\lambda = 15/(x_1 - x_2)$ so that the second equality gives you $y$, i.e. $y = y_2 - 15\frac{y_1 - y_2}{x_1 - x_2}$. So the point of intersection is $(x_2 - 15, y_2 - 15\frac{y_1 - y_2}{x_1 - x_2})$.
You can do the same for $x_1 < x_2$, when you get $x = x_2 + 15$, while for $|x_1 - x_2| < |y_1 - y_2|$ you know that the intersection is either at the top or at the bottom of the square, giving a similar analysis.
This is the problem of finding the intersection of a straight line and a
circle, as commented by J.M.. The more elementary method
using analytical geometry, without rotate or translate the coordinate
axes (which would make the computation easier$^1$), although not being a compact one, is the following (see sketch).
The equation defined by points $
S(N,J)$ and $T(M,I)$ is given by
$$
y-J=m(x-N),\qquad m=\frac{I-J}{M-N}\tag{1}.
$$
The equation of the circle centered at $R$ with radius $d=\overline{RQ}$ is
$$
(x-L)^{2}+(y-H)^{2}=d^{2}.\tag{2}
$$
You need to solve the following system
$$
\left\{
\begin{array}{c}
y-J=m(x-N) \\
(x-L)^{2}+(y-H)^{2}=d^{2},
\end{array}\tag{3}
\right.
$$
which is equivalent to
$$
\left\{
\begin{array}{c}
x=\frac{y-J+mN}{m} \\
\left(\frac{y-J+mN}{m}-L\right)^{2}+(y-H)^{2}=d^{2}.\tag{4}
\end{array}
\right.
$$
Solving the quadratic equation yields (with the help of SWP):
$$
y=\frac{1}{ m^{2}+1 }\left( -mN+Lm+J+m^{2}H\pm \sqrt{\Delta}\right), \tag{5}
$$
where the discriminant is
$$\begin{eqnarray*}
\Delta &=&A+B, \\ \text{with }
A
&=&-m^{4}N^{2}+m^{4}d^{2}-m^{4}L^{2}-m^{2}J^{2}-m^{2}H^{2}+d^{2}m^{2}-2m^{3}NH,
\\
B &=&2Lm^{3}H+2Jm^{2}H+2m^{4}NL-2m^{3}JL+2m^{3}JN.\tag{6}
\end{eqnarray*}$$
The information $\overline{RT}<\overline{RQ}<\overline{RS}$ will define the
signal of the term $ \pm \sqrt{\Delta}$. The coordinates of $Q$ are $O=x,K=y$.
$^1$By making the translation $X=x-L$ and $Y=y-H$, and computing the new coordinates of the points in this $X,Y$ system, the above formulae simplify
(it is equivalent to set $L=H=0$ in them). In the end they should be convert back to the original $x,y$ system.
Best Answer
Here is another way. Choose $t \in [0,1]$, then let $p(t) = (x_1+ t (x_2-x_1), y_1+t(y_2-y_1))$.
Then $p(0) = (x_1,y_1)$, $p(1) = (x_2,y_2)$ and for $t \in (0,1)$, $p(t)$ will be a point in between.
This scheme works even if the two $x$ coordinates are the same.