Temperature distribution for pulse mode Laser operation

I have to repeat the if and elseif statements n times in such a way that for time 1-2 second it should calculate the values (f1,f2,f3,f4) for n=1 and then calculate f1 &f2 for n=2 then it should sum of f1,f2,f3,f4 at n=1 & f1,f2 at n=2 & calculate the equation value for time t-(1-1.5). Then again it should find the values of f1,f2,f3,f4 at n=1 & f1,f2,f3,f4 at n=2 and sum it up & calculate the equation value for time t-(1.5-2). Then again it repeat these values for higher times t-(2-3),Here n=3 & It should calculate all the values in similar fashion mentioned above for n=1,2 & 3 for time t-(2-2.5 & 2.5-3)

clc;close all;

clear all;

K = 51.9;

rho = 7865;

Cp = 472;

T0 = 27;

alpha = K/(rho*Cp);

d = 2.104*10^-3;

a = d/2;

A = (pi/4)*d^2;

DC = 0.5;%(Duty Cycle = 50%)

%v = 50;%(scan speed = 3000mm/min)

v = 1.403*10^-3;

F0 = 200;

I = F0/A;

tint = d/v;%(ti = 1.5sec)

tp = 0.5;%(Laser on,DC = 50%)

tc = 0.5;%(%(Laser off,DC = 50%)

f = 1;%(Frequency in Hz)

z = 0.0001;%(in meters,depth)

t = 0:0.1:3;

for n =1:f

for i = 1:length(t)

if t(i)<=tint

if t(i)<=tp

[delta1(i),f1(i),f2(i)] = nHeating(n,a,z,alpha,tp,tc,t(i));

T(i) = T0+2*(I/K)*(sqrt(pi)/2)*(delta1(i)*(f1(i)-f2(i)));

elseif t(i)>tp && t(i)<=tp+tc

[delta1(i),f1(i),f2(i)] = nHeating(n,a,z,alpha,tp,tc,t(i));

[delta2(i),f3(i),f4(i)] = nCooling(n,a,z,alpha,tp,tc,t(i));

T(i) = T0+2*(I/K)*(sqrt(pi)/2)*((delta1(i)*(f1(i)-f2(i)))-((delta2(i))*(f3(i)-f4(i))));

else t(i)>tint-tp && t(i)<=tint

[delta1(i),f1(i),f2(i)] = nHeating(n,a,z,alpha,tp,tc,t(i));

[delta2(i),f3(i),f4(i)] = nCooling(n,a,z,alpha,tp,tc,t(i));

[delta3(i),f5(i),f6(i)] = nplusoneHeating(n,a,z,alpha,tp,tc,t(i));

T(i) = T0+2*(I/K)*(sqrt(pi)/2)*((delta1(i)*(f1(i)-f2(i)))-(delta2(i)*(f3(i)-f4(i)))+(delta3(i)*(f5(i)-f6(i))));

end

else

[delta3(i),f5(i),f6(i)] = nplusoneHeating(n,a,z,alpha,tp,tc,t(i));

[delta4(i),f7(i),f8(i)] = nplusoneCooling(n,a,z,alpha,tp,tc,t(i));

T(i) = T0+2*(I/K)*(sqrt(pi)/2)*((delta3(i)*(f5(i)-f6(i)))-(delta4(i)*(f7(i)-f8(i))));

end

end

end

T = T(:);

size(T);

plot(t,T,'r')

xlabel('Time in sec');

ylabel('Temperature in Deg.Centigrade');

grid on

function [delta1,f1,f2] = nHeating(n,a,z,alpha,tp,tc,t)

delta1 = sqrt(alpha*(t-((n-1)*(tp+tc))));

f1 = erfc(z/(2*delta1));

f2 = erfc(sqrt(z^2+a^2)/(2*delta1));

end

function [delta2,f3,f4] = nCooling(n,a,z,alpha,tp,tc,t)

delta2 = sqrt(alpha*(t-((n*tp)+((n-1)*tc))));

f3 = erfc(z/(2*delta2));

f4 = erfc(sqrt(z^2+a^2)/(2*delta2));

end

function [delta3,f5,f6] = nplusoneHeating(n,a,z,alpha,tp,tc,t)

delta3 = sqrt(alpha*(t-(n*(tp+tc))));

f5 = erfc(z/(2*delta3));

f6 = erfc(sqrt(z^2+a^2)/(2*delta3));

end

function [delta4,f7,f8] = nplusoneCooling(n,a,z,alpha,tp,tc,t)

delta4 = sqrt(alpha*(t-((n+1)*tp+n*tc)));

f7 = erfc(z/(2*delta4));

f8 = erfc(sqrt(z^2+a^2)/(2*delta4));

end