Dear Ivan,
Thanks for your response. I really appreciate it. Yes, you are right. It has something to do with the usage of lists. So I changed my equation into:
sol = NDSolve[{q'[t] == (a q[t]^b newefr[t] - q[t]) + (a q[t]^b (newefr[t] - newefr1[t])), q[1] == 0.001}, q, {t, 1, 456}]
where:
newefr = Interpolation[newefrain, InterpolationOrder -> 1]; where: newefrain = Table[rainfun[t]*ff, {t, 1, 456}]
similarly to newefr1:
newefr1 = Interpolation[newefrain1, InterpolationOrder -> 1]; where: newefrain1 = Table[rainfun[t]*ff1, {t, 1, 456}]
And it works! Thanks for your help!
Best Regards, Intan