discrete digital filtering in python
Obtaining the system equations

Using Kirchhoff's voltage law:
$$
\begin{align}
  V(t) &amp;= V_{R_1}(t) + V_{R_2}(t) + V_{R_3}(t) + V_I(t) \tag{1}\\
  V_{R_1}(t) &amp;= V_{C_1}(t) \tag{2} \\
  V_{R_2}(t) &amp;= V_{C_2}(t) \tag{3}
\end{align}
$$
And similarly using Kirchhof's current law:
\begin{align}
  I_{R_3}(t) &amp;= I_{R_1}(t) + I_{C_1}(t) \tag{4}\\
             &amp;= I_{R_2}(t) + I_{C_2}(t) \tag{5}
\end{align}

Now using the Voltage-Current relationships of the resistances and capacitors (assuming ideal components):

\begin{align}
  V_{R_1}(t) &amp;= R_1 I_{R_1}(t) \tag{6}\\
  V_{R_2}(t) &amp;= R_2 I_{R_2}(t) \tag{7}\\
  V_{R_3}(t) &amp;= R_3 I_{R_3}(t) \tag{8}\\
  I_{C_1}(t) &amp;= C_1 \frac{d\,V_{C_1}(t)}{dt} \tag{9}\\
  I_{C_2}(t) &amp;= C_2 \frac{d\,V_{C_2}(t)}{dt} \tag{10}
\end{align}

Based on your comments, a current sink with a resistive internal resistance can be assumed. For analysis, it is common to replace such a load with its Norton equivalent circuit (see Norton theorem) which consists of putting the internal resistance (say $R_4$) accross an ideal current sink $I(t)$.
\begin{align}
  V_I(t) &amp;= R_4 I_{R_4}(t) \tag{11} \\
         &amp;= R_4 (I_{R_3}(t) - I_I(t)) \tag{12}
\end{align}

Numerical solution

To discretize the problem to discrete time step $t_n = n \Delta t = n/f_s$, we can make approximations for $(9)$ and $(10)$:
\begin{align*}
  I_{C_1}[n] 
    &amp;\approx C_1 \frac{V_{C_1}[n+1] - V_{C_1}[n]}{\Delta t} \\
    &amp;\approx C_1 f_s \left(V_{R_1}[n+1] - V_{R_1}[n]\right) \\
  I_{C_2}[n] &amp;\approx C_2 \frac{V_{C_2}[n+1] - V_{C_2}[n]}{\Delta t} \\
    &amp;\approx C_1 f_s \left(V_{R_1}[n+1] - V_{R_1}[n]\right)
\end{align*}
This gives us the recurrence equations:
\begin{align*}
  V_{R_1}[n+1] 
    &amp;\approx V_{R_1}[n] + \frac{1}{C_1 f_s} I_{C_1}[n] \\
    &amp;\approx V_{R_1}[n] + \frac{1}{C_1 f_s} \left(I_{R_3}[n] - I_{R_1}\right) \\
    &amp;\approx V_{R_1}[n] + \frac{1}{C_1 f_s} \left(\frac{V_{R_3}[n]}{R_3} - \frac{V_{R_1}[n]}{R_1}\right) \\
  V_{R_2}[n+1] 
    &amp;\approx V_{R_2}[n] + \frac{1}{C_2 f_s} I_{C_2}[n] \\
    &amp;\approx V_{R_1}[n] + \frac{1}{C_2 f_s} \left(I_{R_3}[n] - I_{R_2}\right) \\
    &amp;\approx V_{R_1}[n] + \frac{1}{C_2 f_s} \left(\frac{V_{R_3}[n]}{R_3} - \frac{V_{R_2}[n]}{R_2}\right) 
\end{align*}
And $V_{R_3}[n]$ begin obtained by substituting $(8)$ and $(12)$ into $(1)$:
\begin{align*}
  V_{R_3}[n] 
    &amp;= \frac{V[n] - V_{R_1}[n] - V_{R_2}[n] + R_4 I[n]}{1+R_4/R_3}
\end{align*}

which can readily be implemented in python (with $V(t)$ and $I_I(t)$ being constant for illustration, but the extension to arbitrary function is trivial) with:

# numerical solution
T    = 1.0/fs
Tmax = 20
t = np.arange(0,Tmax,T)

vr3 = np.zeros(len(t))
vr1 = 0
vr2 = 0
for i in np.arange(len(t)-1):
    vr3[i] = (V-vr1-vr2+R4*I)/(1+R4/R3)
    i3  = vr3[i]/R3
    ic1 = i3 - vr1/R1
    ic2 = i3 - vr2/R2
    vr1 = vr1 + ic1/(C1*fs)
    vr2 = vr2 + ic2/(C2*fs)
vr3[-1] = (V-vr1-vr2+R4*I)/(1+R4/R3)


Steady-state solution

As a validation of the above numerical solution, we can check whether it converges to a steady state which can be computed analytically fairly easily (assuming a steady state can be reached, which typically requires $V(t)$ and $I_I(t)$ to be near constant). To do this we note that in the steady state derivative terms go to zero. As a result, $I_{C_1}(t) = I_{C_2}(t) = 0$ and consequently, $I_{R_1}(t)=I_{R_2}(t)=I_{R_3}(t)$. Thus $(1)$ can be written as:
\begin{align*}
  V(t) &amp;= R_1 I_{R_1}(t) + R_2 I_{R_1}(t) + R_3 I_{R_1}(t) + R_4\left(I_{R_1}(t) - I_I(t)\right) \\
     &amp;= \left(R_1+R_2+R_3+R_4\right) I_{R_1}(t) - R_4 I_I(t)
\end{align*}
From which we can solve for $V_{R_3}(t)$:
\begin{align*}
  V_{R_3}(t) &amp;= \frac{R_3}{R_1+R_2+R_3+R_4} \left(V(t) + R_4 I_I(t)\right)
\end{align*}

For illustration, the graph below shows an example for the following parameters:

# parameters
R1 = 2200.0
C1 = 0.001
R2 = 1000.0
C2 = 0.0047
R3 = 5600.0
R4 = 100000.0
V  = 5
I  = 1e-5
fs = 10




As a final note, for voltage source $V(t)$ and current sink $I_I(t)$ which are constant or are from a set of fairly simple functions, the solution could also be obtained analytically with the use of Laplace transforms.