The short answer is that the step response is the inverse Laplace Transform of $\frac{1}{s}G(s)$.
Here are some hints that should help in being able to solve the problem. Since it is likely a homework problem I wanted to provide more details as relevant background and steps toward a solution rather than detail the mathematics for the specific solution. A plot of the resulting step response is included at the end to validate the solution.
First, a review of the simple case of real negative poles on the s-plane along with their inverse (unilateral) Laplace Transform.
As shown, a pole given by the transfer function $H(s) = \frac{1}{s+\alpha}$ has an inverse Laplace transform $x(t) = e^{-\alpha t}$ for $t\ge 0$. This is easy to prove but is also well tabulated in all tables of the Laplace Transforms such as at the bottom of this page: https://en.wikipedia.org/wiki/Laplace_transform
Specifically if $H(s)$ describes the transfer function of the system (output divided by input when the units are $s$), $h(t)$ is the impulse response of the system. As an aside demonstrating one utility of the Laplace Transform, if we restrict $s$ to just be the imaginary axis, we get the frequency response $H(j\omega)$ which in this case would be a low pass filter. Whether $s$ is real or imaginary, we typically refer to the s-domain as the "frequency domain" and all points on the surface are given as pole locations by the transforms of $e^{st}$. So with $h(t)$ being the impulse response, what we see as graphed is the response in time at the output if an impulse was presented at the input. Notice that the pole closest to the imaginary axis takes the longest to decay in its impulse response. This is a key take-away and is the reason we call the negative poles closest to the axis "dominant poles". Like an unwanted guest, they take way too long to go away. Poles further away from the imaginary axis decay much more rapidly in time.
Now consider the summation of these two transfer functions:
$$G(s) = H_1(s) + H_2(s)$$
A summation of two systems in the frequency domain is a parallel process in their signal flow implementation as depicted in the figure below. In contrast the function $H_1(s)H_2(s)$ would be implemented as the cascade of the two systems in series.
Now for the system above represented by the summation of our two systems, the total output due to an impulse at the input would be the sum of the impulse response from each system. Here we see how the pole that was further away from the imaginary axis quickly decays so does not significantly contribute to the final settling time of the output. This settling time (and the overall shape of the response if the poles are far enough apart) is dominated by the pole closest to the imaginary axis, our dominant pole.
Further this combined system is given as:
$$G(s) = \frac{1}{s+b} + \frac{1}{s+a} = \frac{2s+a+b}{(s+a)(s+b)}$$
Notice how adding them together introduced a zero (at $z_1=-(a+b)/2$) but did not change the pole locations. Importantly this shows how given the general combined expression on the right side of the above equation, the process of Partial Fraction Expansions can be used to put this in the form of the parallel (summed) system as described by the left side of the above equation. This is then very easy to solve for the inverse transform, providing the impulse response, here given as:
$$g(t) = e^{-bt} + e^{-at}$$
Step Response
The above shows how to solve for the impulse response, which is the inverse Laplace Transform of the transfer function $G(s)$. To get the step response note that the Laplace Transform of an integrator is $1/s$ (also in the Wikipedia link above and a very important takeaway). Similarly the Laplace Transform of a differentiator is $s$ (less the initial condition). Yes commit these to memory as they will come up again and again and again. So imagine if we prefaced our system under test with an integrator: The integratiion of an impulse is a step! So for any given system, if we simply multiply it's transfer function by $1/s$ (which means putting an integrator in cascade or series with the system), the output defined by the inverse Laplace Transform of that result will be the step response! It's that simple. Taking that further if we multiplied by $1/s^2$ we would get a ramp response, etc.
Solving this again is aided by Partial Fraction Expansion:
$$\frac{1}{s}H_1(s) = \frac{1}{s(s+b)}= \frac{1}{b}\bigg(\frac{1}{s} -\frac{1}{s+b} \bigg)$$
From which we can easily derive the inverse Laplace Transform providing the step response $s(t)$:
$$s(t) = \frac{1}{b}\mathscr{L}^{-1}\bigg(\frac{1}{s} -\frac{1}{s+b} \bigg) = \frac{1}{b}(1-e^{-bt})$$
To easily see how the consideration of the dominant pole still applies, consider that a step at the input to the system is an infinite series of impulses each at time $t$ (and to be accurate each is weighted by $dt$ but that is just a scaling). Each impulse produces it's related impulse response, with the responses from llater impulses (accumulating) with the tails of the responses from previous impulses. The accumulation of the weighted impulse responses is described mathematically as convolution but results in the step response at the output. Therefore the dominant pole would also dominate the resulting step response.
With that background in mind, here are the main points specific to solving the OP's problem:
The step response can reasonably estimated from the dominant pole alone.
For an accurate and complete solution, use Partial Fraction Expansion to factor decompose $G(s)$ into three parallel systems $G_1(s)$, $G_2(s)$, $G_3(s)$, which can then each be easily solved for their inverse Laplace Transforms (this would provide the impulse response of $G(s)$.
The step response for each can be solved as described above by first multiplying by $1/s$ and then solving for the inverse Laplace Transform. The result for this is plotted below with $G_1(s)$ containing the dominant pole $s=-0.1$, $G_2(s)$ containing the weaker pole at $s=-100$ and $G_3(s)$ containing the two repeated poles at $s=-100$ (consistent with Partial Fraction Expansion for repeated poles).
It can be confirmed how closely $G_1(s)$ follows the step response for the system with only the dominant pole given by $1/(s+0.1)$. Also notice the significant difference in time and amplitude scales between the three plots. $G_1(s)$ takes over 50 seconds to settle while $G_2(s)$ and $G_3(s)$ are nearly settled at 50 milliseconds. The sum of all three of these will closely match $G_1(s)$ alone. Further the final values for $G_2(s)$ and $G_3(s)$ are significantly smaller than the final value for $G_1(s)$, but the sum of all three final values will match the final value for combined system $G(s)$.
Comparing the overall step response to that given by the dominant pole highlights the significance of recognizing the dominant pole and how that can be used to simplify most practical problems (I have time now to write this all out from all the time I saved before by recognizing the dominant pole!):
The more significant difference between the two is in the very beginning, so zooming in to show that detail:
The MATLAB/Octave with the control toolbox to plot the step response for $G(s)$ is:
num = -[1 -3]*10000/3;
den = conv([1 0.1], [1 200 10000]);
G = tf(num, den)
step(G);
