Do you mean something more like this?
N = 100; % Total number of particles
one = 1;
alpha = 1;
K = 7; % 7 states
n0 = 4*ones(K,1);
tspan = 0:0.1:5;
[t, n] = ode45(@(t, n) ODEfunc(t, n, K, one, N, alpha), tspan, n0);
plot(t, n),grid
xlabel('time'), ylabel('n')
legend('n1','n2','n3','n4','n5','n6','n7')
function dndt = ODEfunc(~, n, K, one, N, alpha)
S = n(1);
for i=2:K
S = i*n(i) + S;
end
k_bar = S/N;
dndt = zeros(K,1);
dndt(1,:) = alpha*n(1)*k_bar - one*n(1); % Note that this is for state 1.
for i= 2:K-1
dndt(i,:) = alpha*n(i)*(k_bar - i) + one*n(i-1) - one*n(i); % For state 2 to 6
end
dndt(K,:) = alpha*n(K,1)*(k_bar - K) + one*n(K-1,1); % For state 7
end
