How to find the output of a system with the unit-step response and the input?
This can indeed be done. Write the input function $x(t)$ as

$$x(t)=x(-\infty)+\int_{-\infty}^t x^{\prime}(\tau)d\tau=
x(-\infty)+\int_{-\infty}^{\infty} x^{\prime}(\tau)u(t-\tau)d\tau\tag{1}=\\
=x(-\infty)+(x^{\prime}*u)(t)$$

where $x^{\prime}(t)$ is the derivative of $x(t)$, and $u(t)$ is the unit step function. From (1), the system's response $y(t)$ to the input $x(t)$ can be written in terms of the step response $a(t)$:

$$y(t)=a(\infty)x(-\infty)+\int_{-\infty}^{\infty} x^{\prime}(\tau)a(t-\tau)d\tau\tag{2}$$

So, apart from the first constant term, the output signal $y(t)$ is given by the convolution of the derivative of the input signal with the system's step response.

Equation (2) is valid because the system's response to the constant $x(-\infty)$ is $H(0)x(-\infty)$, where $H(0)$ is the system's frequency response at DC:

$$H(0)=\int_{-\infty}^{\infty}h(t)dt=a(\infty)$$

The integral on the right-hand side of (2) is the convolution

$$((x^{\prime}*u)*h)(t)=(x^{\prime}*(u*h))(t)=(x^{\prime}*a)(t)$$

EDIT: The discrete-time case is completely analogous. Write the input sequence as

$$x[n]=x[-\infty]+\sum_{k=-\infty}^n(x[k]-x[k-1])=x[-\infty]+\sum_{k=-\infty}^{\infty}(x[k]-x[k-1])u[n-k]$$

This results in the following expression for the output sequence in terms of the step response:

$$y[n]=x[-\infty]a[\infty]+\sum_{k=-\infty}^{\infty}(x[k]-x[k-1])a[n-k]$$

Again, apart from the constant first term, the output sequence is obtained by convolving the first order difference of the input sequence with the system's step response.