NT = 750;
x1(1,:) = [0, 1];   % coordinate wise would by x1 and y1
x2(1,:) = [0, -1];  % and x2 and y2
%Probably would be possible to use only one x variable with 4 coordinates
%but the solution would be too obfuscated then..
K1 = 1;
K2 = [-1, 1:4, -3];
j = 1:2     %for getting x and y coordinates from the array
k = [-1, 1] %for sign in the x1 and x2 equation
dt = 0.01;  %chosen timestep

for kk = 1:size(K2, 2)
    for ii = 2:NT
        %first index corresponds to x diff, second to y diff in dx
        dx(j) = x1(ii - 1, j) - x2(ii - 1, j);
        r2 = sum(dx .* dx);
        F(j) = dx(flipdim(j, 2)) / r2;
        x1(ii, j) =  x1(ii - 1, j) + K2(kk) * F(j) .* k(j) * dt;
        x2(ii, j) =  x2(ii - 1, j) + K1 * F(j) .* flipdim(k(j), 2) * dt;
    end
    figure('Name', strcat('(K1,K2)=[1,',  int2str(K2(kk)), ']'))
    plot(x1(:,1), x1(:,2), '.', x2(:,1), x2(:,2), 'O');
end

%I didn't include ODE solutions here. It is quite easy anyway. Only a 
% function needs to be defined which output is an array of the solutions
% for y(1) to y(4) and the function fn could be called out smth like that:
%   tt = [0:.01:750];
%   [T,Y] = ode45(@fn,tt,[0 1 0 -1]);
% and that function:
% function dy = fn(y,t)
% dy = zeroes(4,1)
% dy(1) = -K2*(1/r2)*(y(4) - y(2));
% dy(2) = K2*(1/r2)*(y(1) - y(3));
% and etc..