Hi,

(I mailed this answer twice, the first one seems to have not worked properly...)

I believe you made some mistakes in the last 4 integrations. You replace tspan by a range of y (or x) values. What you meant is that you want the solution to stop once it leaves the domain y=[0,4]. However you have an ODE in time, so you cannot know a priori how long it will take for the solution to depart from the saddle point. Here, your y values gets interpreted as a time range (here 4 seconds). When the initial point is very close to the saddle point, the solution departs from it very slowly, so 4 seconds is not sufficient to build the trajectory you are seeking.

If you instead ask for a longer integration time (say 5000 s), you get what you want:

To let the simulations stops when the solution leaves the domain of interest, you should use event functions (see the doc of solve_ivp).

Here is the last part of you code with my corrections:

tspan = np.array([0,5e3])
epsilon = 1e-5
 #TODO: ideally, make event functions so that the simulations are stopped when they leave the domain x=[0,4], y=[0,4] to spare computational time
ic = [x4-epsilon,y4+epsilon]
solution = solve_ivp(f, t_span=tspan, y0=ic, method='RK45', t_eval=None, dense_output=True, args=(a,b,c,d,))
assert solution.success
plt.plot(solution.y[0], solution.y[1],'g--', linewidth=4, zorder=80)
 
ic = [x4+epsilon,y4-epsilon]
solution = solve_ivp(f, t_span=tspan, y0=ic, method='RK45', t_eval=None, dense_output=True, args=(a,b,c,d,))
assert solution.success
plt.plot(solution.y[0], solution.y[1],'m--', linewidth=4, zorder=81)

ic = [x4-epsilon,y4-epsilon]
solution = solve_ivp(f, t_span=-tspan, y0=ic, method='RK45', t_eval=None, dense_output=True, args=(a,b,c,d,))#, rtol=1e-10)
assert solution.success
plt.plot(solution.y[0], solution.y[1], 'y', linewidth=4, zorder=82)
 
ic  = [x4+epsilon,y4+epsilon]
solution = solve_ivp(f, t_span=-tspan, y0=ic, method='RK45', t_eval=None, dense_output=True, args=(a,b,c,d,))#, rtol=1e-10, atol=1E-12)
# assert solution.success
# this last simulation fails because the point reached is too far from the domain of interest and the solution diverges
# however, we already have a nice trajectory inside the required domain
plt.plot(solution.y[0], solution.y[1], 'b', linewidth=4, zorder=83)
# note that the function g=-f is not needed, simply integrate f in negative time instead

Regards,

Laurent