Question: I am trying to fit the parameters of an ODE to data with two dimensions, which should generally be possible, according to the example Fitting multidimensional datasets. This is my failed attempt so far import symfit as sf import numpy as np # data x = np.arange(0,19) data = 10e-4 * np.array([8,10,12,11,10,15,25,37,46,40,43,35,27,14,8,10,13,9,10]) data2 = 10e-3 * np.array([0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0]) # model t, mp, inp = sf.variables('t, mp, inp') tau, w = sf.parameters('tau, w') model = {sf.D(mp, t): (-mp + inp * w) / tau } # fitting ode_model = sf.ODEModel(model, initial={inp: data2[0], t: 0.0, mp: data[0]}) fit = sf.Fit(ode_model, inp=data2, t=x, mp=data) fit_result = fit.execute()  It seems I'm not correctly defining inp in the model definition. I'm getting the following error TypeError: got an unexpected keyword argument 'inp' . I suspect that I'm making a mistake in providing inp as named_data to sf.Fit() as ist does not appear to be an independent variable in the model documentation sf.Fit() This is the full error message: /usr/local/lib/python3.10/dist-packages/symfit/core/fit.py in __init__(self, model, objective, minimizer, constraints, absolute_sigma, *ordered_data, **named_data) 372 # Bind as much as possible the provided arguments. 373 signature = self._make_signature() --> 374 bound_arguments = signature.bind_partial(*ordered_data, **named_data) 375 376 # Select objective function to use. Has to be done before calling /usr/lib/python3.10/inspect.py in bind_partial(self, *args, **kwargs) 3191 Raises `TypeError` if the passed arguments can not be bound. 3192 """ -> 3193 return self._bind(args, kwargs, partial=True) 3194 3195 def __reduce__(self): /usr/lib/python3.10/inspect.py in _bind(self, args, kwargs, partial) 3173 arguments[kwargs_param.name] = kwargs 3174 else: -> 3175 raise TypeError( 3176 'got an unexpected keyword argument {arg!r}'.format( 3177 arg=next(iter(kwargs)))) TypeError: got an unexpected keyword argument 'inp'  Could someone help? Thank you so so much :-)

Answer: inp has to be defined as an expression and integrated into the expression d mp / dt. To do so, the data for inp has to be fit so as to reproduce data2. Since data2 looks like a square wave, a Fourier series is used to fit the data. The following is the code implementing these ideas: import symfit as sf import numpy as np from functools import reduce import matplotlib.pyplot as plt from scipy.optimize import curve_fit from scipy.integrate import solve_ivp def square_wave_symbolic(x, L, shift, n=8): return reduce(lambda a, b: a + b, [sf.sin(2*(1+2 *k) *sf.pi * (x+shift) / L)/(1+2* k) for k in range(n)]) def square_wave_numeric(x, L, shift, n=8): return reduce(lambda a, b: a + b, [np.sin(2*(1+2 *k) *np.pi * (x+shift) / L)/(1+2* k) for k in range(n)]) # data x = np.arange(0,19) data = 10e-4 * np.array([8,10,12,11,10,15,25,37,46,40,43,35,27,14,8,10,13,9,10]) data2 = 10e-3 * np.array([0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0]) # Fit square wave alpha0 = 5E-3 beta0 = 5E-3 L0 = 18 shift0 = -4 n = 9 plt.figure() plt.title('Approximation of inp using Fourier Series') result = curve_fit(lambda x, alpha, beta, shift, L: alpha + beta * square_wave_numeric(x, L, shift, n=n), x, data2, p0=(alpha0, beta0, shift0, L0), full_output=True) alpha, beta, shift, Lf = result[0] plt.scatter(x, data2, label='Original data') plt.plot(x, alpha + beta * square_wave_numeric(x, Lf, shift, n=n), label='Fit data') plt.legend() # model t, mp = sf.variables('t, mp') tau, w = sf.parameters('tau, w') inp = alpha + beta * square_wave_symbolic(t, Lf, shift, n=n) model = {sf.D(mp, t): (-mp + inp * w) / tau } # fitting ode_model = sf.ODEModel(model, {t: 0.0, mp: data[0]}) fit = sf.Fit(ode_model, t=x, mp=data) fit_result = fit.execute() print(fit_result.params) w = fit_result.params['w'] tau = fit_result.params['tau'] # Verify parameters by solving numerically def f(t, x, w, tau): mp = x[0] inp = alpha + beta * square_wave_numeric(t, Lf, shift, n=n) return np.array([(-mp + inp * w) / tau]) t0 = x[0] tf = x[-1] t_eval = np.linspace(t0, tf, 200) ode_result = solve_ivp(f, (t0, tf), (data[0],), t_eval=t_eval, args=(w, tau), method='Radau') plt.figure() plt.scatter(x, data, label='Original data') plt.plot(ode_result.t, ode_result.y[0], label='ODE solver data') plt.legend()  I end up with the following value for fit_result.params: OrderedDict([('tau', 1.191491195628205), ('w', 3.4675371620653133)])  The following are the plot for inp variable fit:  ODE fit plot:  Notes:  The square wave function describing inp is not smooth. I am not sure about where the data comes from, but if this is just some test data, the real data may have to be fit using a different function.  It would make sense to do the ODE and inp function fit together. Doing it using a minimize function is feasible. I am not sure how this could be done with symfit as I am not familiar with this library.