Defining 2N ODEs with 2N variables when using ode45

I am new to Matlab and am using it to complete a university project. I am trying to solve a system of 2N ODEs for the position of N vortices in a 2D fluid. They interact with each other and have a circulation strenth Γ, and each vortex has a location given by the and coordinate for each index (See below the system of ODEs). I am able to set up the equations explicitly with a couple of 'for' loops but am struggling to figure out how to define each equation when it is in a 'function' environment and without a nested 'for' loop.

Each vortex has an initial position for x,y which I have each stored in 2 Nx1 vectors and I want to define each differential equation in terms of x and y as shown in the screenshot below (note the prime on the sum denotes an ommision of the α= β term). I also am trying to use ode45 to solve these and then plot the output which I am able to do, it's just the defining of the ODEs where I am stuck.

Hope this makes sense and thanks in advance for any help.

[James Tursa]

Show us the loops you have written. It should be a simple matter of putting these loops into a function that is suitable for ode45( ).

[William]

function dzdt = vortex_pos(t,z)
%% Point Vortices
N = 4;
z = zeros(N,2); %first column x, second column y
x = z(:,1);
y = z(:,2);
z = [x;y];

dzdt = zeros(2*N,1); 
x_sum = zeros(N,1);
y_sum = zeros(N,1);

%% ODEs
for i = 1:N
    for j = setdiff(1:N, i)
        x_sum(i) = x_sum(i) + (-1/(2*pi))*(gamma(j)*(y(i)-y(j))/((x(i)-x(j))^2 + (y(i)-y(j))^2));
        y_sum(i) = y_sum(i) + (1/(2*pi))*(gamma(j)*(x(i)-x(j))/((x(i)-x(j))^2 + (y(i)-y(j))^2));
    end
    dzdt(i) = x_sum(i);
    dzdt(N+i) = y_sum(i);
end

end

This is my attempt at writing the function for the ODEs. I just use z to store both the x and y variables for later use with ode45. I also added x_sum because for some reason I couldn't use the dzdt term in the loop in that place. The code does run when I use ode45 but I just get NaN for all elements of z.

[Torsten]

You overwrite the painfully computed vector of solution variables z by zeros:

function dzdt = vortex_pos(t,z)
%% Point Vortices
N = 4;
z = zeros(N,2); %first column x, second column y

So all the l_ij in the denominators of your summands become zero and the expression NaN.

[William]

Thanks I have sorted it now. Also what do you mean by painfully? Do you mean by way of the for loops?

[Torsten]

I meant that ODE45 spent so much effort computing the solution and you just carelessly overwrite it :-)

function dzdt = vortex_pos(t,z)
  N = numel(z)/2;
  x = z(1:N);
  y = z(N+1:end);
  dzdt = zeros(2*N,1);
  %% ODEs
  for i = 1:N
    for j = setdiff(1:N, i)
        dzdt(i) = dzdt(i) + (-1/(2*pi))*(gamma(j)*(y(i)-y(j))/((x(i)-x(j))^2 + (y(i)-y(j))^2));
        dzdt(N+i) = dzdt(N+i) + (1/(2*pi))*(gamma(j)*(x(i)-x(j))/((x(i)-x(j))^2 + (y(i)-y(j))^2));
    end
  end
end