%runtime: 26 seconds dt=0.01; rcore=0.1; dz0=0; Klist=[1 -1 1 -1]; nvortex=length(Klist); zlist=[-1-0.5i -1+0.5i -0.5-0.5i -0.5+0.5i]; i=nvortex+1; for x=-1:0.1:1; for y=-1:0.1:1; zlist(i)=x+1i*y; i=i+1; end; end; n=length(zlist); for i=1:1:n; for j=1:1:10; paths(i,j)=zlist(i); end; end; for t=1:1:85; for i=1:1:n; dz(i)=0; for j=1:1:nvortex; if(zlist(i)!=zlist(j)); r2=abs(zlist(i)-zlist(j))^2; dz(i)=dz(i)+1i*Klist(j)*(zlist(j)-zlist(i))/r2*(1-exp(-r2/rcore^2)); end; end; end; zlist=zlist+dt*(1.5*dz-0.5*dz0); dz0=dz; for i=1:1:n; temp=cat(2,zlist(i),paths(i,:)); temp(11)=[]; paths(i,:)=temp; end; if(t==1); hold off; else; hold on; end; for i=1:1:n; plot(real(paths(i,:)),imag(paths(i,:))); end; axis([-1 1 -1 1]); axis square; hold off; end;