See below for a demonstration, i had to assume some parameters. However the plotting princaple remains the same.
% these parameters were missing from the OP's question, hence assumed to be 1
alpha = 1;
qs = 1;
K = 1;
Ti = 1;
% parameters from the OP's question
y = 0.01;
t = y:10*y:1000*y;
x = 0:0.001:1;
% initialize matrix T
T = zeros(length(x),length(t));
for l = 1:length(x)
for k = 1:length(t)
a = 60*sqrt(4*alpha*t(k));
b = a/sqrt(pi);
c = qs/K;
d = exp(-(x(l)./a).^2);
T(l,k) = c*(b*d-x(l).*erfc(x(l)./a))+Ti;
end
end
ts = t(:,4);
% make the full plot
figure
plot(x,T)
grid on
%% only plot some selected 't' values
% first pick a few t values
T_pick = [1 5 10 20 50 75 100];
% in the plot command, use the 2nd dimension to plot the selected values
figure
plot(x, T(:,T_pick))
grid on
