### Find intersections numerically -- Using NumPy

Intersects of 2 functions Assuming function $f$ and $g$, scipy.optimize.fsolve would give an elegant solution: from scipy.optimize import fsolve def f(xy): x, y = xy return (y - x ** 2, y - np.sin(x)) fsolve(f, [1.0, 1.0]) gives array([0.87672622, 0.76864886]) is the intersect of $y=x^2$ and $y=\sin(x)$. Intersects of 2 arrays There is a function called numpy.intersect1d that returns the sorted, uniq values in both arrays, and reduce(np.intersect1d, list_of_arrays) gives intersect of all arrays in the list. Index can also be returned by return_indices=True option. If the 2 arrays are data values, using np.flatnonzero(np.abs(arr1-arr2))<tol) seems to be a good solution; this method returns the indices of values in arr1 and arr2 which are closed by tol . Another method is based on idea that there is a change in sign of $\lim_{x\to x_0^+}f(x_0)-g(x_0)$ and $\lim_{x\to x_0^-}f(x_0)-g(x_0)$, therefore, if $\mathrm{d}x$ is small enough, this method would give a very good e