You can also write it out more explicitly. The equation reads:
[ 1,pt2 ; pt3,pt4 ] * [ pt5 ; 1 ] = [ pt5 + pt2 ; pt3.*pt5 + pt4 ]
Since each of those terms is a scalar, you can compute them for all t at the same time using element-wise multiplication:
t = linspace(0,100);
pt2 = t+(1/2)*exp(-3*t)-(1/2)*exp(-t);
pt3 = (3/2)*(exp(-t)-exp(-3*t));
pt4 = 1-(3/2)*exp(-3*t)+(1/2)*exp(-t);
pt5 = (t-4)/3;
y_mat = (1./t) .* [ pt5 + pt2 ; pt3.*pt5 + pt4 ];
plot(t,y_mat)
This might be a bit more verbose, but I don't think it's any less readable than other solutions. And it is much more efficient: 0.0571 ms, versus 483.3 ms (syms solution) and 0.681 ms (loop solution), for a t with 500 elements.
(Note that multiplying by 1./t uses implicit singleton expansion. This works in MATLAB R2016b and newer. For older versions of MATLAB, use bsxfun.)