loop to repeat steps and plot on same figure

Below is the code displayed in full. Is there a way to complete all the same processes but in a shorter script using foor loops to repeat the plots and combine where there is repetition of the same equations? Hopefully this makes sense :)

Qr = 20.83;    % m3/s
B = 13;        % m 
h = 1.3;       % m 
V = Qr/(B*h);  % m/s
XS_Area=B.*h;     % m2

%% 2 instantaneous releases 
t1= 5*60;         % s 
t2= 10*60;        % s
Mo= (5*10)/2;     % kg
M1=Mo./XS_Area;   % kg/m2

C1=(M1./sqrt(4.*3.141.*E.*t1) ).*exp(-(X1-V.*t1).^2./(4.*E.*t1));
C2=(M1./sqrt(4.*3.141.*E.*t2) ).*exp(-(X1-V.*t2).^2./(4.*E.*t2));
C3=C1+C2;

plot(X1,C3,'g-'),xlabel('X(m)'),ylabel('C_o(kg/m^3)'),grid minor,xlim([0 1500]);
hold on;

%% 5 instantaneous releases;
t01=2*60;
t02=4*60;
t03=6*60;
t04=8*60;
t05=10*60;

Mo5=(5*10)/5;       % kg
M15=Mo5./XS_Area;   % kg/m2

C01=(M15./sqrt(4.*3.141.*E.*t01) ).*exp(-(X1-V.*t01).^2./(4.*E.*t01));
C02=(M15./sqrt(4.*3.141.*E.*t02) ).*exp(-(X1-V.*t02).^2./(4.*E.*t02));
C03=(M15./sqrt(4.*3.141.*E.*t03) ).*exp(-(X1-V.*t03).^2./(4.*E.*t03));
C04=(M15./sqrt(4.*3.141.*E.*t04) ).*exp(-(X1-V.*t04).^2./(4.*E.*t04));
C05=(M15./sqrt(4.*3.141.*E.*t05) ).*exp(-(X1-V.*t05).^2./(4.*E.*t05));

C5=C01+C02+C03+C04+C05;
plot(X1,C5,'b-'),xlabel('X(m)'),ylabel('C_o(kg/m^3)'),grid minor,xlim([0 1500]);
grid minor;
hold on;

%% 10 instantaneous releases; 
t001=1*60;
t002=2*60;
t003=3*60;
t004=4*60;
t005=5*60;
t006=6*60; 
t007=7*60;
t008=8*60;
t009=9*60;
t010=10*60;

Mo10=(5*10)/10;        % kg
M1_10=Mo10./XS_Area;   % kg/m2

C001=(M1_10./sqrt(4.*3.141.*E.*t001) ).*exp(-(X1-V.*t001).^2./(4.*E.*t001)); 
C002=(M1_10./sqrt(4.*3.141.*E.*t002) ).*exp(-(X1-V.*t002).^2./(4.*E.*t002));
C003=(M1_10./sqrt(4.*3.141.*E.*t003) ).*exp(-(X1-V.*t003).^2./(4.*E.*t003));
C004=(M1_10./sqrt(4.*3.141.*E.*t004) ).*exp(-(X1-V.*t004).^2./(4.*E.*t004));
C005=(M1_10./sqrt(4.*3.141.*E.*t005) ).*exp(-(X1-V.*t005).^2./(4.*E.*t005));
C006=(M1_10./sqrt(4.*3.141.*E.*t006) ).*exp(-(X1-V.*t006).^2./(4.*E.*t006));
C007=(M1_10./sqrt(4.*3.141.*E.*t007) ).*exp(-(X1-V.*t007).^2./(4.*E.*t007));
C008=(M1_10./sqrt(4.*3.141.*E.*t008) ).*exp(-(X1-V.*t008).^2./(4.*E.*t008));
C009=(M1_10./sqrt(4.*3.141.*E.*t009) ).*exp(-(X1-V.*t009).^2./(4.*E.*t009));
C010=(M1_10./sqrt(4.*3.141.*E.*t010) ).*exp(-(X1-V.*t010).^2./(4.*E.*t010));

C10 = C001+C002+C003+C004+C005+C006+C007+C008+C009+C010; 
plot(X1,C10,'k-'),xlabel('X(m)'),ylabel('C_o(kg/m^3)'),grid minor,xlim([0 1500]);
grid minor;
hold off;

[Rik]

Numbered variables should be replaced by indexing of an array.

See this discussion.

Answer 1

[David Hill]

%write simple script to run function
Qr = 20.83;    
B = 13;       
h = 1.3; 
E = ?;%unknown
X1=(1:.1:10)';%whatever
T={60*(5:5:10),60*(2:2:10),60*(1:10)};
figure;hold on;
for k=1:length(T)
    yourFunction(Qr,B,X1,E,h,T{k});
end


%X1 needs to be column vector, t needs to be row vector, assume E is scalar
function yourFunction(Qr, B, X1, E, h, t)
  V = Qr/(B*h);
  XS_Area=B.*h;
  M=(5*10)/length(t)/XS_Area;
  C=sum(M./sqrt(4*pi*E*t).*exp(-(X1-V*t).^2./(4*E*t)),2);
  plot(X1,C);
end

[Ben Bawtree]

Thank you! that has done it :)