c=3e8;M=5;P=10 % since you did not define them I gave P and M these values
N=[20:20:640]
gamma_array=[]
xj_array=[]
big_dzdx=zeros(101,32) % 32 out of the 640/20 you want N to sweep through
for j=1:1:numel(N)
i=1;
for x=0:c/N:c
if x>=(P*c)
dzdx(:,i)=(M*c/(1-P)^2)*(2*P/c-2*x/c^2);
else
dzdx(:,i)=(M*c/P^2)*(2*P/c-2*x/c^2);
end
i=i+1;
end
x=[0:c/N:c];
big_dzdx(:,j)=dzdx
%%Using vortex_lattice to get gamma and xj
[gamma, xj] = vortex_lattice(dzdx, x, alpha, vinf, rhoinf, c);)
gamma_array(:,j)=[gamma_array; gamma]
xj_array=[xj_array; xj]
end
