clear all
close all
clc

% network size
N = 50;

% global coupling parameter values used in manuscript.
%WE = 1.5;
%WE = 2.05;
WE = 2.115;
%WE = 2.25;
for WE = [1.5,2.05,2.115,2.25]
% other parameters 
theta = 1;
p = 0.2;
a = 5;
tau1 = 1;
tau2 = 5;
T = 1*10^3; 

%% create a weak coupling/diagonally dominant matrix 
epsilon = 1e-3;
WEE = eye(N) + epsilon*rand(N,N);
for j = 1:N
WEE(j,:)=WEE(j,:)*WE/sum(WEE(j,:));    
end
    
%% 
options = odeset('AbsTol',1e-14,'RelTol',1e-14);
int = rand(3*N,1);
tic
[t,y] = ode45(@(t,y) WC(WEE,tau1,tau2,p,a,theta,N,y),0:0.1:T,int,options); 
[t,y] = ode45(@(t,y) WC(WEE,tau1,tau2,p,a,theta,N,y),0:0.1:T,y(end,:)',options); 
toc
[tsn2,ysn2]= ode45(@(t,y) WC(WE,tau1,tau2,p,a,theta,1,y),0:0.1:T,[y(1,1),y(1,N+1),y(1,2*N+1)]',options); 
save(sprintf('WEAK_COUPLING_sim_WE=%3.3f.mat',WE));


%% 
figure
title(sprintf('WE=%3.3f.mat',WE))
subplot(1,2,1)
plot(y(:,1:10),y(:,N+1:N+10))
xlabel('E')
ylabel('I')
subplot(1,2,2)
plot(t,y(:,1:N))
xlabel('t')
ylabel('E(t)')
xlim([0,1000])
end