Question: In this minimal example I want to solve the basic integrator ODE dT/dt = k[F(t) - T(t)] where F(t) is a forcing vector, which is selected to be a square wave. I have implemented two procedures to solve the ODE: Using Scipy's solve_ivp and using Gekko. The former takes 6 milliseconds to run, the latter takes 8 seconds to run. Am I doing something wrong with Gekko? Are there some additional parameters that I can tune to improve performance? I have omitted the testing part of the code, but I have compared both solutions to the analytic solution, and they are both accurate within a 3-4 significant digits. import numpy as np from scipy.integrate import solve_ivp from gekko import GEKKO def step_problem_forcing(time_arr, param): T_low = param['T_low'] T_high = param['T_high'] t_break_l = param['t_break_l'] t_break_r = param['t_break_r'] forcing = np.full(len(time_arr), T_low) idxs_pulse = np.logical_and(time_arr >= t_break_l, time_arr <= t_break_r) forcing[idxs_pulse] = T_high return forcing def step_problem_generate(tmin, tmax): tAvg = (tmin+tmax) / 2 tDur = tmax - tmin return { 'k' : 10 ** np.random.uniform(-3, 3), 'T_low' : np.random.uniform(20, 50), 'T_high' : np.random.uniform(20, 50), 'T0' : np.random.uniform(20, 50), 't_break_l' : np.random.uniform(tmin, tAvg - 0.1*tDur), 't_break_r' : np.random.uniform(tAvg + 0.1*tDur, tmax) } def ode_scipy_step_solver(time_arr, param: dict): f_this = lambda t: np.interp(t, time_arr, param['forcing']) def _dxdt(t, x, param: tuple): # if (t < tmin) or (t > tmax): # raise ValueError(f"{t} is out of bounds [{tmin}, {tmax}]") k, = param return k*(f_this(t) - x) k = param['k'] sol = solve_ivp(_dxdt, (time_arr[0], time_arr[-1]), [param['T0'], ], t_eval=time_arr, args=([k, ],), rtol=1.0E-6, method='RK45') return sol['y'].T def ode_gekko_step_solver(time_arr, param: dict) -> np.ndarray: m = GEKKO() # create GEKKO model # Define variables m.time = time_arr T = m.Var(value=param['T0']) F = m.Param(value=param['forcing']) k = m.Const(value=param['k']) # t = m.Param(value=m.time) # equations m.Equation(T.dt() == k * (F - T)) m.options.IMODE = 7 # dynamic simulation m.solve(disp=False) # solve locally (remote=False) return T.value tmin = 0 tmax = 10 time_arr = np.linspace(tmin, tmax, 1000) param = step_problem_generate(tmin, tmax) param['forcing'] = step_problem_forcing(time_arr, param) # Takes 6.8ms to run rezScipy = ode_scipy_step_solver(time_arr, param) # Takes 8.5s to run rezGekko = ode_gekko_step_solver(time_arr, param) 

Answer: Switch to remote=False to solve locally and avoid the overhead due to network communication to the public server. Switching to IMODE=4 solves all 1000 time steps simultaneously instead of sequentially and can improve solve time. Here are 10 trials with the 1 original and 9 modified versions. I'm on Ubuntu Linux. Timing will likely be different based on the computer. # IMODE=7, Original Script ODE Scipy Time: 0.014239072799682617 Gekko Time: 8.92684006690979 # IMODE=4, Remote=False ODE Scipy Time: 0.2553427219390869 Gekko Time: 0.12177252769470215 ODE Scipy Time: 0.0017511844635009766 Gekko Time: 0.11760997772216797 ODE Scipy Time: 0.0032286643981933594 Gekko Time: 0.11528849601745605 ODE Scipy Time: 0.06327152252197266 Gekko Time: 0.12140154838562012 ODE Scipy Time: 0.01037144660949707 Gekko Time: 0.12077593803405762 ODE Scipy Time: 0.0036330223083496094 Gekko Time: 0.11592268943786621 ODE Scipy Time: 0.1100163459777832 Gekko Time: 0.11649894714355469 ODE Scipy Time: 0.0032815933227539062 Gekko Time: 0.1154778003692627 ODE Scipy Time: 0.0015087127685546875 Gekko Time: 0.1158149242401123  Here's the modified script: import numpy as np from scipy.integrate import solve_ivp from gekko import GEKKO import time def step_problem_forcing(time_arr, param): T_low = param['T_low'] T_high = param['T_high'] t_break_l = param['t_break_l'] t_break_r = param['t_break_r'] forcing = np.full(len(time_arr), T_low) idxs_pulse = np.logical_and(time_arr >= t_break_l, time_arr <= t_break_r) forcing[idxs_pulse] = T_high return forcing def step_problem_generate(tmin, tmax): tAvg = (tmin+tmax) / 2 tDur = tmax - tmin return { 'k' : 10 ** np.random.uniform(-3, 3), 'T_low' : np.random.uniform(20, 50), 'T_high' : np.random.uniform(20, 50), 'T0' : np.random.uniform(20, 50), 't_break_l' : np.random.uniform(tmin, tAvg - 0.1*tDur), 't_break_r' : np.random.uniform(tAvg + 0.1*tDur, tmax) } def ode_scipy_step_solver(time_arr, param: dict): f_this = lambda t: np.interp(t, time_arr, param['forcing']) def _dxdt(t, x, param: tuple): # if (t < tmin) or (t > tmax): # raise ValueError(f"{t} is out of bounds [{tmin}, {tmax}]") k, = param return k*(f_this(t) - x) k = param['k'] sol = solve_ivp(_dxdt, (time_arr[0], time_arr[-1]), [param['T0'], ], t_eval=time_arr, args=([k, ],), rtol=1.0E-6, method='RK45') return sol['y'].T def ode_gekko_step_solver(time_arr, param: dict) -> np.ndarray: m = GEKKO(remote=False) # create GEKKO model # Define variables m.time = time_arr T = m.Var(value=param['T0']) F = m.Param(value=param['forcing']) k = m.Const(value=param['k']) # t = m.Param(value=m.time) # equations m.Equation(T.dt() == k * (F - T)) m.options.IMODE = 4 # dynamic simulation m.solve(disp=False) # solve locally (remote=False) return T.value tmin = 0 tmax = 10 time_arr = np.linspace(tmin, tmax, 1000) param = step_problem_generate(tmin, tmax) param['forcing'] = step_problem_forcing(time_arr, param) st = time.time() rezScipy = ode_scipy_step_solver(time_arr, param) print(f'ODE Scipy Time: {time.time()-st}') st = time.time() rezGekko = ode_gekko_step_solver(time_arr, param) print(f'Gekko Time: {time.time()-st}')  Scipy solve_ivp is a specialized ODE solver while Gekko solves higher index DAEs, mixed integer optimization, regression, and machine learning models. It uses a simultaneous method and only has basic adaptive step lengths to control error. I recommend solve_ivp if it is an ODE problem.