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:

1. 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.

2. 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.