For part (1), since you know $x_1 = x_3$ and $x_2 = x_4$, and compressing or stretching the lines horizontally won't change the $y$-coordinate of the intersection, you can just solve the problem for $x_1 = x_3 = 0$ and $x_2 = x_4 = 1$ in full generality. The exact $x$-coordinates are actually irrelevant.
More formally, define new coordinates $x', y'$, where $y = y'$ and $x = x_1 + (x_2 - x_1) x'$. Thus we can write $x'_1 = x'_3 = 0$ and $x'_2 = x'_4 = 1$. The equations of your lines in the new coordinates (dropping the redundant prime marks on $y$) are $$y - y_1 = (y_2 - y_1) x'$$ and $$y - y_3 = (y_4 - y_3) x'.$$
Rearrange the second equation to $$x' = \frac{y - y_3}{y_4 - y_3}$$ and substitute into the first equation to get \begin{align*}y - y_1 &= \frac{(y_2 - y_1) (y - y_3)}{y_4 - y_3} \\ \left( 1 - \frac{y_2 - y_1}{y_4 - y_3} \right)y &= y_1 - \frac{(y_2 - y_1) y_3}{y_4 - y_3} \\
(y_1 - y_2 - y_3 + y_4) y &= (y_4 - y_3) y_1 - (y_2 - y_1) y_3 \\
y &= \frac{y_1 y_4 - y_2 y_3}{y_1 - y_2 - y_3 + y_4}. \end{align*}
Two sanity checks: $y$ is invariant under interchange of $(y_1, y_2)$ with $(y_3, y_4)$, and in the case $y_1 = y_4 = 0$ and $y_2 = y_3 = 1$ we get $y = \frac{1}{2}$, as we would expect.
For part (2), I don't think you could improve anything. If you knew the $x$-coordinates, then you could use quadratic or cubic splines, but with unknown $x$-coordinates, this is impossible.
Best Answer
You can just check the source code.
Leading to minpack. Hybrd and hybrj are essentially the same, but hybrd uses forward differences to compute the jacobian whereas hybrj requires the user to provide the jacobian. They use Powell's method, with the modifications described in the previous link to minpack.