**Parametric solution**

If the sequences {x1,y1} and {x2,y2} define arbitrary (x,y) curves, rather than y(x) curves, we need a parametric approach to finding intersections. Since it's not entirely obvious how to do so, and because @unutbu's solution uses a defunct interpolator in SciPy, I thought it might be useful to revisit this question.

```
import numpy as np
from numpy.linalg import norm
from scipy.optimize import fsolve
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
x1_array = np.array([1,2,3,4,6,8,9])
y1_array = np.array([4,12,7,1,6.3,8.5,12])
x2_array = np.array([1.4,2.1,3,5.9,8,9,23])
y2_array = np.array([2.3,3.1,1,3.9,8,9,11])
s1_array = np.linspace(0,1,num=len(x1_array))
s2_array = np.linspace(0,1,num=len(x2_array))
# Arguments given to interp1d:
# - extrapolate: to make sure we don't get a fatal value error when fsolve searches
# beyond the bounds of [0,1]
# - copy: use refs to the arrays
# - assume_sorted: because s_array ('x') increases monotonically across [0,1]
kwargs_ = dict(fill_value='extrapolate', copy=False, assume_sorted=True)
x1_interp = interp1d(s1_array,x1_array, **kwargs_)
y1_interp = interp1d(s1_array,y1_array, **kwargs_)
x2_interp = interp1d(s2_array,x2_array, **kwargs_)
y2_interp = interp1d(s2_array,y2_array, **kwargs_)
xydiff_lambda = lambda s12: (np.abs(x1_interp(s12[0])-x2_interp(s12[1])),
np.abs(y1_interp(s12[0])-y2_interp(s12[1])))
s12_intercept, _, ier, mesg \
= fsolve(xydiff_lambda, [0.5, 0.3], full_output=True)
xy1_intercept = x1_interp(s12_intercept[0]),y1_interp(s12_intercept[0])
xy2_intercept = x2_interp(s12_intercept[1]),y2_interp(s12_intercept[1])
plt.plot(x1_interp(s1_array),y1_interp(s1_array),'b.', ls='-', label='x1 data')
plt.plot(x2_interp(s2_array),y2_interp(s2_array),'r.', ls='-', label='x2 data')
if s12_intercept[0]>0 and s12_intercept[0]<1:
plt.plot(*xy1_intercept,'bo', ms=12, label='x1 intercept')
plt.plot(*xy2_intercept,'ro', ms=8, label='x2 intercept')
plt.legend()
print('intercept @ s1={}, s2={}\n'.format(s12_intercept[0],s12_intercept[1]),
'intercept @ xy1={}\n'.format(np.array(xy1_intercept)),
'intercept @ xy2={}\n'.format(np.array(xy2_intercept)),
'fsolve apparent success? {}: "{}"\n'.format(ier==1,mesg,),
'is intercept really good? {}\n'.format(s12_intercept[0]>=0 and s12_intercept[0]<=1
and s12_intercept[1]>=0 and s12_intercept[1]<=1
and np.isclose(0,norm(xydiff_lambda(s12_intercept)))) )
```

which returns, for this particular choice of initial guess [0.5,0.3]:

```
intercept @ s1=0.4761904761904762, s2=0.3825944170771757
intercept @ xy1=[3.85714286 1.85714286]
intercept @ xy2=[3.85714286 1.85714286]
fsolve apparent success? True: "The solution converged."
is intercept really good? True
```

This method only finds one intersection: we would need to iterate over several initial guesses (as @unutbu's code does), check their veracity, and eliminate duplicates using `np.close`

. Note that `fsolve`

may falsely indicate successful detection of an intersection in the return value `ier`

, which is why the extra checking is done here.

Here's the plot for this solution: