if you treat the diodes as having memoryless non-linear volt-amp characteristics and treat the capacitor as linear and having memory, you can use Euler's backward method to represent the capacitors and everything else are static Kirchoff equations (with some nonlinearity, you need to represent the back-to-back diodes accurately - the standard diode equation might not be good enough).
okay, first let's deal with the parallel diodes. a single diode volt-amp characteristic is
$$ i_D(t) = I_D \left( e^{\frac{q \ v_D(t)}{k T}} -1 \right) $$
$v_D(t)>0$ means forward bias of the diode and a significant current $i_D(t)>0$ flows. when $v_D(t)<0$, it is reversed biased and the most current that will flow is $i_D(t)=-I_D<0$ . $I_D$ is the reverse saturation current and is very tiny.
$q$ is the elementary charge (sometimes called the "electron charge", but is a positive number), $k$ is the Boltzmann constant, and $T$ (in this section) is the absolute temperature of the PN junction of the diode (about 293K). $\frac{kT}{q}$ is a voltage about $\tfrac{1}{40}$ volt. when $v_D(t)=\frac{kT}{q}$ the forward diode current is still approximately $+I_D$, which is extremely small. the exponential curve doesn't visibly kick in until $v_D(t) \approx 25 \cdot \frac{kT}{q} = 0.65$ volt.
keeping the polarities defined in the same direction but reversing the diode results in
$$ i_D(t) = -I_D \left(e^{\frac{-q \ v_D(t)}{k T}} -1\right) $$
putting the two in parallel (assuming the same $I_D$), the volt-amp description is
$$\begin{align}
i_D(t) & = I_D \left(e^{\frac{q \ v_D(t)}{k T}} -1\right) - I_D \left(e^{\frac{-q \ v_D(t)}{k T}} -1\right) \\
\\
& = I_D \left(e^{\frac{q \ v_D(t)}{k T}} - e^{\frac{-q \ v_D(t)}{k T}}\right) \\
\\
& = 2 I_D \sinh\left(\frac{q \ v_D(t)}{k T}\right) \\
\end{align} $$
inverting the volt-amp equation:
$$\begin{align}
v_D(t) &= \frac{kT}{q}\operatorname{arcsinh}\left(\frac{i_D(t)}{2 I_D} \right) \\
&= \frac{kT}{q}\ln\left(\frac{i_D(t)+\sqrt{i_D(t)^2 + (2I_D)^2}}{2 I_D} \right) \\
\end{align} $$
probably the diodes have some contact resistance, and there is a fudge factor $\eta \approx 1$ or $2$ scaling the voltage, which will make the equation look like
$$ v_D(t) = \frac{\eta kT}{q}\operatorname{arcsinh}\left(\frac{i_D(t)}{2 I_D} \right) + R_D \, i_D(t) $$
now, probably, rather than trying to find (from measured curves) the values of $I_D$ and $R_D$ and, the $\eta\frac{kT}{q}$ factor, you will probably just have a power series to represent the diodes. since the equation has odd-symmetry, your power series will have only odd-order terms.
$$ i_D(t) = a_1 v_D(t) + a_3 (v_D(t))^3 + a_5 (v_D(t))^5 + \dots $$
you might not even need the inner terms and might be able to get away with just two terms, like
$$ i_D(t) = a_1 v_D(t) + \dots a_9 (v_D(t))^9 $$
ignoring terms in between. this will make your computations cheaper.
that's all i want to say about the diodes. you will have some memoryless model for the diode pair:
$$ i_D(t) = f\big( \, v_D(t) \,\big) $$
or the inverse function
$$ v_D(t) = g\big( \, i_D(t) \,\big) $$
that you will need to nail down a decent approximation of.
Euler's backward method:
the capacitor has fundamental volt-amp characteristics:
$$\begin{align}
i_C(t) &= C \frac{d \, v_C(t)}{dt} \\
&= C \lim_{\Delta t \to 0} \frac{v_C(t) - v_C(t-\Delta t)}{\Delta t} \\
\\
&\approx C \frac{v_C(t) - v_C(t-\Delta t)}{\Delta t} \qquad \text{for very small } \Delta t \\
\end{align} $$
or looking at it from the perspective of capacitor voltage
$$\begin{align}
v_C(t) &= \frac{1}{C} \int\limits_{-\infty}^{t} i_C(u) \, du \\
&= \frac{1}{C} \left( \int\limits_{-\infty}^{t-\Delta t} i_C(u) \, du + \int\limits_{t-\Delta t}^{t} i_C(u) \, du \right) \\
&= v_C(t-\Delta t) + \frac{1}{C} \int\limits_{t-\Delta t}^{t} i_C(u) \, du \\
&\approx v_C(t-\Delta t) + \frac{1}{C} \int\limits_{t-\Delta t}^{t} i_C(t) \, du \\
\\
&= v_C(t-\Delta t) + \frac{1}{C} i_C(t) \, \Delta t \qquad \text{for very small } \Delta t \\
\end{align} $$
both approximations say exactly the same thing and demonstrate the meaning we have when we differentiate between memoryless devices and those with memory. the ideal diode is memoryless because, in the ideal, the relationship between the voltage and current at its terminals depends only on their values at the present. the ideal capacitor is a device with memory because the relationship between the present voltage and current at its terminals depends on a voltage the capacitor remembers from the very recent past.
so let the tiny "differential" time $\Delta t$ be our sampling period or $\Delta t = \frac{1}{f_\text{s}}$. then we need evaluate $t$ at only integer sampling periods:
$$ t = n \cdot \Delta t $$
then
$$\begin{align}
v_C(t) &= v_C(t-\Delta t) + \frac{\Delta t}{C} i_C(t) \\
\\
v_C(n \Delta t) &= v_C(n \Delta t-\Delta t) + \frac{\Delta t}{C} i_C(n \Delta t) \\
\\
&= v_C((n-1) \Delta t) + \frac{\Delta t}{C} i_C(n \Delta t) \\
\\
v_C[n] &= v_C[n-1] + \frac{\Delta t}{C} i_C[n] \\
\end{align}$$
at a specific time $t$, for every capacitor, you can represent that capacitor as a (briefly) constant voltage source (the previous sample's voltage) in series with a resistor having value $R_C=\tfrac{\Delta t}{C}$. and we can express that voltage source in series with a resistor as a Thevenin source having its own Norton equivalent whenever it is convenient to express it as such. (it is convenient to do that for the 470 pF cap in the feedback path.)
let $C_1$=100nF, $R_1$=1K, $C_2$=470pF, $R_2$=10K, assuming the op-amp is ideal, the load resistance is not salient. the op-amp in negative feedback maintains virtual equality across its input terminals. the current flowing through $R_1$ into $C_1$ is
$$\begin{align}
i_1(t) &= \frac{1}{R_1}(v_\text{in}(t) - v_{C1}(t)) \\
&= \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}(t) - v_{C1}(t-\Delta t)) \\
\\
i_1[n] &= \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}[n] - v_{C1}[n-1]) \\
\end{align}$$
now the feedback current must equal $i_1(t)$. here, you want to express the source/resistance model of $C_2$ as a Norton equivalent. then the two diodes, the 10K resistor $R_2$, and the Norton resistance ($\tfrac{\Delta t}{C_2}$ are in parallel with the Norton current source of $v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t}$) must be placed in parallel and an aggregate volt-amp characteristic must be known. you want voltage as a function of current.
$$ v_D(t) = g\big( \, i_D(t) \, \big) $$
that is the volt-amp characteristic of two diodes in parallel with two known resistances. so the feedback current will be divided into two currents, the known Norton equivalent current from the known $C_2$ voltage of the previous sample period and whatever flows into your non-linear device which is two diodes in parallel with two resistances.
$$ i_1(t) = -v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} + i_D(t) $$
or
$$\begin{align}
i_D(t) &= i_1(t) + v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} \\
&= \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}(t) - v_{C1}(t-\Delta t)) + v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} \\
\end{align}$$
$$\begin{align}
v_D(t) &= g\big( \, i_D(t) \, \big) \\
&= g\big( i_1(t) + v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} \big) \\
&= g\bigg( \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}(t) - v_{C1}(t-\Delta t)) + v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} \bigg) \\
\end{align}$$
in discrete time, it's
$$v_D[n] = g\bigg( \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}[n] - v_{C1}[n-1]) + v_{C2}[n-1]\tfrac{C_2}{\Delta t} \bigg) $$
and the output is
$$\begin{align}
v_\text{out}[n] &= v_D[n] + v_\text{in}[n] \\
&= g\bigg( i_1[n] + v_{C2}[n-1]\tfrac{C_2}{\Delta t} \bigg) + v_\text{in}[n] \\
&= g\bigg( \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}[n] - v_{C1}[n-1]) + v_{C2}[n-1]\tfrac{C_2}{\Delta t} \bigg) + v_\text{in}[n] \\
\end{align}$$
and you must update your two capacitor states for the next sampling period:
$$\begin{align}
v_{C1}[n] &= v_{C1}[n-1] + \frac{\Delta t}{C_1} i_1[n] \\
&= v_{C1}[n-1] + \frac{\Delta t}{C_1} \left( \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}[n] - v_{C1}[n-1]) \right) \\
\\
v_{C2}[n] &= v_D[n] \\
\end{align}$$
