dt=0.01; 
rcore=0.1; 
dz0=0; 
Klist=[1,-1,1,-1];
zlist=[-1-0.5i,-1+0.5i,-0.5-0.5i,-0.5+0.5i]; 
i=length(Klist)+1;

for x=-5:0.1:5;
    for y=-5:0.1:5; 
        zlist(i)=x+1i*y; 
        i=i+1;
    end;
end;

n=length(zlist);

for i=1:n;
    for j=1:10;
        paths(i,j)=zlist(i);
    end;
end;

for t=1:85;
    for i=1:n; dz(i)=0;
        for j=1:length(Klist);
            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;
    hold on;
    for i=1:n; 
        temp=cat(2,zlist(i),paths(i,:)); 
        temp(11)=[]; paths(i,:)=temp;
        plot(real(temp),imag(temp));
    end;
    axis([-1,1,-1,1]);
    axis square;
    hold off;
    drawnow
end;