If I recall correctly, back in R2015a, fplot did not yet apply to symbolic expressions.
That symsum() to infinity does not find a closed form solution so it leaves it as a symsum. If you use matlabFunction() on the result, then it leaves it in terms of symsum(), which is not a good thing as symsum is not defined for numeric parameters. You can attempt to fplot() the resulting handle, but it takes too long for any practical purposes. If you substitute in specific t values and run the symsum with those, it takes too long for any practical purposes.
I think for practical purposes you are going to need to truncate the infinite sums.
You also have the problem that you have factorial(k*mu) but mu = 1.5 so for odd integers k*mu is non-integral which is a problem for factorial. You have to switch to gamma.
TERMS = 30;
syms k t
I_0=1; mu=1.5; lambda=0.3; Gamma=0.1;
t=0:1:20;
inner1 = (-lambda *t.^(mu)).^k./(2.^((k*(k-1)*mu)/2).*gamma(k*mu+1));
inner2 = ((-lambda).^(k-1).*(t.^(k*mu)))./(2.^((k*(k-1)*mu)/2).*gamma(k*mu+1));
inner1sum = symsum(inner1, k, [0, TERMS]);
inner2sum = symsum(inner2, k, [1, TERMS]); %notice different bounds
y = I_0 * inner1sum + Gamma * inner2sum;
Y = double(y);
plot(t, Y)
grid on
xlabel('time t');
ylabel('I(t)');
