How to find the intersection time of a parameterized curve with a shape?

I have a curve parameterized by time that intersects a shape (in this case just a rectangle). Following this elegant suggestion, I used shapely to determine where the objects intersect, however from there on, I struggle to find a good solution for when that happens. Currently, I am approximating the time awkwardly by finding the point of the curve that is closest (in space) to the intersection, and then using its time stamp.

But I believe there should be a better solution e.g. by solving the polynomial equation, maybe using the root method of a numpy polynomial. I'm just not sure how to do this, because I guess you would need to somehow introduce tolerances as it is likely that the curve will never assume exactly the same intersection coordinates as determined by shapely.

Here is my code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Ellipse
from matplotlib.collections import LineCollection
from shapely.geometry import LineString, Polygon


# the parameterized curve
coeffs = np.array([
    [-2.65053088e-05, 2.76890591e-05],
    [-5.70681576e-02, -2.69415587e-01],
    [7.92564148e+02, 6.88557419e+02],
])
t_fit = np.linspace(-2400, 3600, 1000)
x_fit = np.polyval(coeffs[:, 0], t_fit)
y_fit = np.polyval(coeffs[:, 1], t_fit)
curve = LineString(np.column_stack((x_fit, y_fit)))

# the shape it intersects
area = {'x': [700, 1000], 'y': [1300, 1400]}
area_shape = Polygon([
    (area['x'][0], area['y'][0]),
    (area['x'][1], area['y'][0]),
    (area['x'][1], area['y'][1]),
    (area['x'][0], area['y'][1]),
])

# attempt at finding the time of intersection
intersection = curve.intersection(area_shape).coords[-1]
distances = np.hypot(x_fit-intersection[0], y_fit-intersection[1])
idx = np.where(distances == min(distances))
fit_intersection = x_fit[idx][0], y_fit[idx][0]
t_intersection = t_fit[idx]
print(t_intersection)

# code for visualization
fig, ax = plt.subplots(figsize=(5, 5))
ax.margins(0.4, 0.2)
ax.invert_yaxis()

area_artist = Rectangle(
    (area['x'][0], area['y'][0]),
    width=area['x'][1] - area['x'][0],
    height=area['y'][1] - area['y'][0],
    edgecolor='gray', facecolor='none'
)
ax.add_artist(area_artist)

points = np.array([x_fit, y_fit]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
z = np.linspace(0, 1, points.shape[0])
norm = plt.Normalize(z.min(), z.max())
lc = LineCollection(
    segments, cmap='autumn', norm=norm, alpha=1,
    linewidths=2, picker=8, capstyle='round',
    joinstyle='round'
)
lc.set_array(z)
ax.add_collection(lc)

ax.autoscale_view()
ax.relim()

trans = (ax.transData + ax.transAxes.inverted()).transform
intersection_point = Ellipse(
    xy=trans(fit_intersection), width=0.02, height=0.02, fc='none',
    ec='black', transform=ax.transAxes, zorder=3,
)
ax.add_artist(intersection_point)

plt.show()

And just for the visuals, here is what the problem looks like in a plot:

from How to find the intersection time of a parameterized curve with a shape?