R//C circuit in time domain in Python

I am trying to simulate a resistance in parallel of a capacitor with a varying current as input (current vs time, with a fast changing pace). I want to do this in Python in the time domain using the exact solution of the equation. For validation, I put the differential equation in Simulink and compare the results.

Equation resolution

Equation: dU / dt + U / (C*R) = I/C

Regular solution: dU / dt + U / (CR) = 0 => u(t) = kexp(-t/(RC))+B

Particular solution: I inject u(t) = kexp(-t/(RC))+B in dU / dt + U / (CR) = I/C => B = RI

Boundary condition: u(0) = Ureff giving: k = Ureff - RI

So general solution (work in discharge, charge, rest...): u(t)=(Ureff - RI) * exp(-t/(RC)) + RI

Python

• I determine the indexes where the current changes (if current passes from 5 to 10 at between index 9 and10, I will save the index 10)
• I have a loop going through the time, with the equation: Utotal[i] = (Uref-RI)np.exp(-(time-offset_time) / (RC))+RI
• Uref and offset time are updated using the index where the current change (using a if). Uref takes as value the Utotal[i-1] and offset_time the value at time[i]. I=current[i]

The issue

It works well when the current does not change too fast (whatever R and C) or when R and C are below 1 and the current is changing fast. When I use the same current profile, but with values of R and C above 1, the Utotal[i] totally overshoots (it follows RI).

Here the voltage vs time at R=1 and C=1 (Python in red, Matlab in purple) Here at R=1 and C=10

Question

Is it an issue with my solution of the equation (though it is working in some case) ?

Is it an issue with Uref taken at Utotal[i-1] when a change in current happens at i (and if the current changes too fast and with a high R*C, it does not work anymore) ?

Thank you

itprorh66 is correct, I interpolated all the input using a small sampling rate and it is much more accurate. it takes forever though and some parts are still not good, maybe the algo can be improved. Then I used a formula found in the litterature using discretization: Urc(t) =Urc(t-1)exp(-dt/(RC)) + RI(t)*(1-exp(-dt/(RC)) and it works fast and as Simulink and close to measurements (now remains to understand how they found the formula above, not really clear besides it doing discretization)