Hi @Arina_Belova. Unfortunately, I did not yet have time to look into this issue in more detail, but sometimes the analytical solutions (generated by using method='exact') can contain terms that create a division by zero for some parameter values. A typical example would be when you have two time constants (e.g. tau_1 and tau_2 in your example) with the same values – or of course when one of the time constants is zero. I don’t see anything obvious like this in your equations, though.
As a workaround, you can change the integration method to e.g. method='rk2', this will probably make the error go away.
For debugging, the technique I mention in this comment could be helpful: ZeroDivisionError: float division - #2 by mstimberg