How can I convert a step response that is described by the sum of two exponents to an IIR filter analytically? If a system has the following step response:

$$\begin{align}
y(t) &amp;= 1 - \Big(A e^{−t/T_1} + (1-A) e^{−t/T_2}\Big) \\
 &amp;= 1 - A e^{−t/T_1} - (1-A) e^{−t/T_2} \\
\end{align}$$

and I know the values of $A$, $T_1$ and $T_2$, can I transform this to a discrete-time IIR filter analytically?  If there is no analytical solution, is it possible to say how many poles/zeros the IIR filter should have?

I've done it numerically by using the Prony's method on the impulse response, but I'm hoping there is a more direct way.
You can use the step-invariant transformation, which makes sure that the step response of the discrete-time system matches the step response of the continuous-time system at the sampling instants. I suppose this is what you mean by "exact solution".

Sampling the given continuous-time step response at $t=nT$ gives the desired step response of the discrete-time system:

$$g[n]=y(nT)=\left[1-Ae^{-nT/T_1}-(1-A)e^{-nT/T_2}\right]u[n]\tag{1}$$

where $u[n]$ is the unit step.
The $\mathcal{Z}$-transform of $(1)$ is

$$G(z)=\frac{1}{1-z^{-1}}-\frac{A}{1-e^{-T/T_1}z^{-1}}-\frac{1-A}{1-e^{-T/T_2}z^{-1}}\tag{2}$$

The frequency response of the desired discrete-time system is then given by

$$H(z)=(1-z^{-1})G(z)\tag{3}$$

After performing some basic algebra and using $a=e^{-T/T_1}$ and $b=e^{-T/T_2}$ we obtain

$$H(z)=\frac{[1-b+A(b-a)]z^{-1}+[a(b-1)+A(a-b)]z^{-2}}{1-(a+b)z^{-1}+abz^{-2}}\tag{4}$$

As a result, the IIR filter with transfer function $(4)$ has a step response that equals the given continuous-time step response $y(t)$ at the sample instants $t=nT$: