need help with a program to plot the phase of 100 oscillators as they change in time. The phases of the oscillators are also coupled

I am trying to plot the phase of 100 oscillators versus time, each with a different frequency. The average of the frequency is supposed to be non-zero. The phases of the oscillators are coupled. The oscillators have a periodicity of 2pi.

theta=2*pi*rand(10); %multiplying 2pi by 10x10 array of random values 
%between 0 and 1 to generate 100 random initial values of theta between
%0 and 2pi
omega=10*rand(10); %similarly multiplying 10 by 10x10 array of random 
%values between 0 and 1 to generate 100 random values of omega between 
%0 and 10

if sum(sum(omega))/100 == 0
    omega=10*rand(10); 
else disp(omega) 
end %this conditional statement is meant to assure the average of all omega
%is a non-zero value

theta_circmean = atan(.001*(sum(sum(sin(theta))))/(sum(sum(cos(theta)))));
%this is the circular mean of all theta values

for n = linspace(1,10,10)
    for m=linspace(1,10,10)
        t=linspace(1,100,100);
        theta(n,m)=theta(n,m)+5*sin(theta_circmean-theta(n,m)); %coupling
        %oscillators to eachother using circular mean
        theta_t(n,m,t) = theta(n,m)+omega(n,m)*t; %pairing theta(n,m) with 
        %omega(n,m) to create create 100 functions of theta_t in time
        theta_t(n,m,t)=mod(theta_t(n,m,t)+2*pi,2*pi);%making theta_t periodic
    end
end

plot(t,theta_t)

[Bish Erbas]

What is your question?

[Jessica Johnson]

So what I have above does not give me the plot (time versus phase) that I am looking for. I am sure it is for a few reasons: 1) how do I couple all the theta values properly? This was my attempt:

theta(n,m)=theta(n,m)+5*sin(theta_circmean-theta(n,m))

2) how do I make theta periodic? i defined it and used mod() to try to do so:

theta_t(n,m,t) = theta(n,m)+omega(n,m)*t; 
theta_t(n,m,t)=mod(theta_t(n,m,t)+2*pi,2*pi)

These are the most pressing questions. Thanks.