In a broad solution, you need to reverse the convolution of your system's output.

Assuming it's a linear continuous system, and it's Single Input and Single Output (SISO), you do the Laplace transform of the signal L{y(t)}=Y(s), obtain the Laplace transform of the input L{x(t)}=X(s), and then obtain the transfer function of the system: H(s)=Y(s)/X(s), you must be aware the transfer function of the step is 1/s, therefore: H(s)=Y(s)/(1/s) => H(s)=sY(s), then you do the inverse Laplace transform: L-¹{H(s)}=L-¹{sY(s)}, which, depending on your system, may require partial fraction expansion. By the end you have h(t) (got a bit lazy here since y(t) is not known, but the step function is very well known).

Of course I made a bunch of assumptions about your system, if your system has discrete steps, the Z transform is of interest, with its own caveats mind you. Then there are filters and other numerical approximations for a reverse convolution.