Is there a linear filter whose step response is an S-curve?
Sums of harmonic cosines

The Hann window is a sum of a rectangular function and a truncated cosine. As an impulse response it has the qualities you asked for:

Figure 1. Hann window as impulse response, and the corresponding step response. The step response is valued $x - \frac{\sin(2\pi x)}{2\pi}$ over $0 \le x \le 1$, and is constant outside that range, because the impulse response is zero outside that range. The step response and its first two derivatives are continuous, because the impulse response and its first derivative are continuous.

A filter with Hann impulse response can be constructed as the (normalized) sum of a truncated integrator filter and the real part of the lowest-non-zero-frequency bin of a sliding discrete Fourier transform, both of which can be calculated recursively as truncated IIR filters.

The filters have poles on the z-plane unit circle (they are marginally stable) and require periodical refreshing of their state to cut accumulation of round-off errors. For fixed point input the truncated integrator can be implemented with extended bit depth so that no rounding errors accumulate and refreshing its state is not necessary.

Octave source for the composite filter:

r = 20;
w = 2*pi/r;
b0 = zeros(1, r+1);
b0(1) = 1;
b0(r+1) = -1;
a0 = [1, -1];
b1 = zeros(1, r+2);
b1(1) = -1; b1(2) = cos(w); b1(r+1) = 1; b1(r+2) = -cos(w);
a1 = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
ir = filter(b0, a0, impulse) + filter(b1, a1, impulse);
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/r #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");


Figure 2. Hann impulse response from the first Octave script.

It is possible to get rid of the unnecessary one-sample delay without additional cost by multiplying the complex transfer function of the sliding DFT bin by $e^{-i\omega}$ before taking the real part:

r = 20;
w = 2*pi/r;
b0 = zeros(1, r+1);
b0(1) = 1;
b0(r+1) = -1;
a0 = [1, -1];
b1 = zeros(1, r+2);
b1(1) = -cos(w); b1(2) = 1; b1(r+1) = cos(w); b1(r+2) = -1;
a1 = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
ir = filter(b0, a0, impulse) + filter(b1, a1, impulse);
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/r #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");


Figure 3. Undelayed Hann impulse response from the second Octave script.

There are other window functions with harmonic cosine terms, and filters with them as impulse responses could be implemented similarly.

B-splines

Piece-wise polynomials such as B-splines are another alternative for the impulse response shape. Cascaded integrator–comb (CIC) filters can be used to implement them. A triangle function would be the minimal piece-wise polynomial function (and B-spline) to meet your demands. The implementation is very simple: A truncated integrator in cascade with an identical truncated integrator. This comes from the favorable property of B-spline functions that they are each a convolution of an impulse by rectangular function multiple times. The number of continuous derivatives is incremented by each convolution. The earlier notes about stability apply.

Figure 4. A triangle impulse response and the corresponding step response.

Octave source:

r = 20;
b = zeros(1, r+1);
b(1) = 1;
b(r+1) = -1;
a = [1, -1];
length = 2*r+10;
impulse = eye(1, length);
ir = filter(b, a, filter(b, a, impulse));
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/(r*r) #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");


Figure 4. Triangular impulse response from the Octave script.