Ah I see.
I am not too familiar with finding values at specific time t in this way though.
Do you need the posterior distribution for each of the input1, input2…? Or do you just need the values of input1, input2?
If it is the former, you will need to do a pm.Gamma or pm.Normal for each of the input1, input2.
pm.Gamma(input1, alpha=alpha, beta=beta)
If it is the latter, you will just need to solve_ivp or odeint outside of pm?