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: [image: image.png] 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
