As unclean as it may be, you can place the dopant drive in figure generation in a loop with different end times
for Times = [301,601,901]
%% Dopant Drive in
T_drive= 900; %anneal temp in celsius
Dop = 10.5; %preexponential diffusion term for P in cm2/s
Ea = 4.28e4; %activation energy
D_drive=Dop*exp(-1*Ea/(T_drive+273)); %calculating the diffusivity at specified T
dx = x(2)-x(1);
dt=(0.5*(dx)^2)/D_drive;
x_d=[-1e-7:1e-7:100e-7];
time_d = (0:dt:Times); %time interval for annealing, seconds
I'd recommend putting the legend in after this loop, assuming you know which time you want.
As I say, not the cleanest, but it works.
