Actually, your solution is not that far off if you keep in mind that the indefinite integral is only determined up to the addition of a constant. However, here is my way to compute the integral:
First off, we need to clarify: The Heaviside function I'll use is defined by
$$H(x) = \left\lbrace \begin{array}{cl} 1 & , ~ x \geq 0 \\ 0 & , ~ x < 0 \end{array}\right.$$
An indefinite integral is just a function that can be written as a definite integral with a variable as the upper bound:
$$\left( \int x \cdot H(x-a) ~ dx \right)(x_0) = \int\limits_a^{x_0} x \cdot H(x-a) ~ dx$$
Now we just have to divide the upper term into two cases:
If $x_0 \geq a$ then $H(x-a) = 1$ for all $x \in [a,x_0]$ and therefore the (right) intgral equals
$$\int\limits_a^{x_0} x~dx = \frac{1}{2}(x_0)^2 - \frac{1}{2} a^2 = \frac{1}{2} (x_0 - a)(x_0+a)$$
For $x_0 <a$ we have $H(x-a) = 0$ for all $x \in [x_0, a]$ and therefore the integral equals $0$.
So if we merge both cases one may write
$$\int\limits_a^{x_0} x \cdot H(x-a) ~ dx = \frac{1}{2} (x_0-a)(x_0 + a)H(x_0 - a)$$
which is the solution you obtained using Mathematica.
