How can I convert a step response that is described by the sum of two exponents to an IIR filter analytically?
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$: