The solution is not continuous. This appears cleary in what the OP writes "The correct impulse response should be : y(t) = exp(-t/tau)*Heaviside(t)/tau" (a claim you can find in any textbook) and in your solution too.
From the theoritical solution provided by the OP, one can easily check it doesn't verify y(0)=0 because Heaviside(0) is undefined. Taking the limit (from the right for a causal system) at t=0+, doesn't solve the issue
limit(exp(-t/tau)/tau*Heaviside(t), t=0, right)
Thus y(0)=0, y(t) verifies the ode D(y)(t) + y(t)/tau = Dirac(t)/tau in the open half line (0, +oo), but y(t) is discontinuous at t=0.
Mybe what puzzles you in the solution I got is that absence of the Heaviside function?
I agree it's quite dusturbing, but this is a Maple problem:
laplace(Heaviside(t), t, p);
invlaplace(%, p, t)
To find a solution which verifies both the ode and the initial condition on can use a smoothed ode where the Dirac term is replaced by a function f(t) such that int(f(t), t=0..+oo)=1 and f(t) is almost equal to 0 in the open interval (epsilon..+oo) where epsilon is a small strictly positive number.
use Statistics in f := PDF('Gamma' (lambda, 2), t) end use:
This function has a peak close to 0, which is higher the smaller lambda is (lambda > 0)
Name it "an engineer trick" or "a physicist's tinkering" if you want, or consider this as a kind of "weak solution"