Could I use a function to achieve this task?

Hi guys,

My code here generates a wave from a source (wave a). I want to create a diagonal line of N sources (Starting from wave a, ending at wave n). As you can see from my code, I've manually added a second wave b in diagonally, changing [x,y] manually as well - where it says %offsets secondary source).

I want to know if, instead of adding each source manually, I can create a function to do this that will generate N amounts of sources in a diagonal line, and I want there to be a time delay between each one. Wave a propagates first, then wave b, then wave c, a bit like a domino effect.

clear
%define the total time for simulation
T_total=5;
%set the computation increment
dt=0.05;
%calculate the number of steps
ntstep=T_total/dt+1;
%define Input Waves
alphaa=2;  %velocity of wave a
fa=1;       %frequency of wave a
ampa=1;     %amplitude of wave a

alphab=2;  %velocity of wave b
fb=1;       %frequency of wave b
ampb=1;     %amplitude of wave b

%Calc wave parameters using input
Ta=1/fa;    %period of wave a
wa=2*pi*fa; %angular frequency of wave a
la=Ta*alphaa;%wavelength of wave a
ka=2*pi/la; %wavenumber of wave a

Tb=1/fb;    %period of wave b
wb=2*pi*fb; %angular frequency of wave b
lb=Tb*alphab;%wavelength of wave b
kb=2*pi/lb; %wavenumber of wave b

x=(-20:0.1:20);        %define the x coordinates 
y=(-20:0.1:20);        %define the y coordinates
t=zeros(length(x),length(y));        %start time at each point in x by y grid (zero everywhere)
%setup the figure
figure(1);
%label axes
hold on
[X,Y]=meshgrid(x,y); %create the mesh: matrix X contains x-by-y values of x, Y contains x-by-y values of y
R1=sqrt(X.^2+Y.^2);  %define the distances to all points on the x-by-y grid.
R2=sqrt((X-5).^2+(Y+5).^2); %Offsets secondary source 

%loop over nsamp
for n=1:ntstep
    t=t+dt; %time
    arga=wa*t-ka*R1; %argument of wave a
    a=ampa*sin(arga)./max(1,sqrt(R1));    %amplitude of wave a
    a(R1>t.*alphaa)=0; %causality condition for a
    a(R1<t.*alphaa-la/2)=0; %limit length of input wave a to one half cycle

      argb=wb*t-kb*R2; %argument of wave a
      b=ampb*sin(argb)./max(1,sqrt(R2));    %amplitude of wave a
      b(R2>t.*alphab)=0; %causality condition for a
      b(R2<t.*alphab-lb/2)=0; %limit length of input wave a to one half cycle

      %set up the figure
      figure(1);
      hold off;

      surf(X,Y,a+b); %plot the sum of the waves a and b at all locations defined by [X,Y]
      axis square; %set square axes
      shading INTERP; %interpolate the shading
      caxis([0,1]);  %define the limits of the colour bar
      c=colorbar;    %draw the colour bar
      ylabel(c,'Amplitude (m)'); % label the colour bar
      view (2); %set the view to aerial
      hold on;
      %label axes
      xlabel('Distance x (m)');
      ylabel('Distance y (m)');
      %limit axes
      xlim([x(1), x(length(x))]);
      ylim([x(1), y(length(y))]);
      zlim([0,ampa]);
      title(strcat('time =', num2str(t(1)),' s'));
      %pause (dt); % can pause each increment of the simulation for slower
      %run-time

end