Hi Jeremy,
The equation for the convolution integral in the Question is incorrect. The argument of u inside the integral should be tau, not t. Also, when dealing with the delta function it's better to evaluate the convolution integral over the entire real line to ensure we integrate through the delta funciton, i.e., if the uppler limit is t, then the integrand "stops" at dirac(0), which is problemlatic.
So
syms tau t real
syms u(t) d(t) y(t)
y(t) = int(u(tau)*d(t-tau),tau,-inf,inf,'Hold',true)
Nwo define the signals u(t) and d(t) and evaluate y(t)
u(t) = heaviside(t);
d(t) = dirac(t);
y(t) = (release(subs(y(t))))
That's a funny looking expression, but we can see that y(t<0) = 0, y(0) = 1/2, and y(t>0) = 1, which is heaviside(t) (using the default definiton for the heaviside(0)). Plot it to verify y(t) is the step function.
fplot(y(t))
ylim([-0.5 1.5])
Recall that the general rule for convolution with the dirac function is that the convolution just returns the other signal, as just shown for example, but this can get a little sticky when the other signal has discontinuous jumps, like the step function.
We can try the sympref, as @Walter Roberson did, to assume that heaviside(0) = 1 as typical for the unit step function, but I don't think it will matter. When t = 0, the convolution integral is
y0 = int(u(tau)*d(0-tau),tau,-inf,inf,'Hold',true)
I think the rule, at least in engineering math, for this situation (heaviside discontinuous at t =0) is that the integral evaluates to the average of heaviside(0-) and heaviside(0+), i.e., the values of heaviside(t) just to the left and right of the origin, which will always be 1/2.
