Hi everyone, I have an issue with integrate.solve_ivp. In a matlab example with matlabs function [~,x] = ode45(@(t,x) g(t,x,a,b,c,d),tspan,ic), I am able to create this output, visible the 4 solution lines from the saddle point. [cid:0b1985d9-c2c1-4011-b08b-4c5ef63d70c5] The problematic parts are the last 4 blocks calculating solutions from the saddle point at [1,1], the unstable node [0,0] and the stable nodes [0,c] and [a,0]. This function is defined as follows: def lotkavolterra_variant(a=3, b=2, c=2, d=1, xmin=0, xmax=4, ymin=0, ymax=4, stepsize=.25): ## Lotka-Volterra system with intraspecific competition # Skript Figure 4.27 x = np.arange(xmin, xmax+stepsize, stepsize) y = np.arange(ymin, ymax+stepsize, stepsize) yval = np.arange(ymin, ymax+stepsize*0.1, stepsize*0.1) # evaluate curve 4 times more often to get smoothness [X, Y] = np.meshgrid(x, y) dx = X*(a - X - b*Y) dy = Y*(c - Y - d*X) f = lambda t, x,a,b,c,d : [x[0]*(a - x[0] - b*x[1]), x[1]*(c - x[1] - d*x[0])] g = lambda t, x,a,b,c,d : [-x[0]*(a - x[0] - b*x[1]), -x[1]*(c - x[1] - d*x[0])] plt.figure(figsize = (13,13)) plt.quiver(X, Y, dx, dy, angles = 'xy', color='#0086b3', width=0.0015) plt.grid() plt.xlabel('x') plt.ylabel('y') plt.axis([xmin,xmax,ymin,ymax]) tspan=[0,50] x0stepsize = .0001 for x0 in np.arange(0, .01+x0stepsize, x0stepsize): ic = [x0, .01-x0] solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,)) plt.plot(solution.y[0], solution.y[1],'r') x0stepsize = .1 for x0 in np.arange(0, 10+x0stepsize, x0stepsize): ic = [x0, 10-x0] solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,)) plt.plot(solution.y[0], solution.y[1],'r') x4 = (a - b*c)/(1 - b*d) y4 = (c - a*d)/(1 - b*d) plt.plot( 0, 0 ,'go',fillstyle='none', ms=15, zorder=100)[0].set_clip_on(False) plt.plot (0, c ,'yo',fillstyle='full', ms=15, zorder=100)[0].set_clip_on(False) plt.plot( a, 0 ,'mo',fillstyle='full', ms=15, zorder=100)[0].set_clip_on(False) plt.plot(x4,y4 ,'bo',fillstyle='none', ms=15, zorder=100)[0].set_clip_on(False) epsilon = 0.00001 ic = [x4-epsilon,y4+epsilon] solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,)) plt.plot(solution.y[0], solution.y[1],'g--', linewidth=2, zorder=80) ic = [x4+epsilon,y4-epsilon] solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,)) plt.plot(solution.y[0], solution.y[1],'m--', linewidth=2, zorder=81) ic = [x4-epsilon,y4-epsilon] solution = solve_ivp(g, [ymin, ymax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,)) plt.plot(solution.y[0], solution.y[1], 'y', linewidth=2, zorder=82) ic = [x4+epsilon,y4+epsilon] solution = solve_ivp(g, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,)) plt.plot(solution.y[0], solution.y[1], 'b', linewidth=2, zorder=83) While Matlab is able to compute them, the Scipy implementation failed with lots of variations of methods, first_step, max_step (I even tried to set it to 0.0001), rtol and atol adjustments. The result is always like this, so the computaion stops after about 40 steps and a short distance [cid:cbbec946-971a-4768-965c-086199212999] unless I set epsilon to something quite high like 0.1. Maybe I'm missing something, so I hope someone can help me resolve to a setting, which calculates the solution for those ode's like the Matlab implementation. Thanks in advance and kind regards