clear all
close all
clc



 clear all
close all
clc


% global coupling parameter values used in manuscript.
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; 
%Lattice Connectivity 
Q = 17;
N = Q^2;
W = zeros(Q,Q);

theta = 1;
p = 0.2;
a = 5;
tau1 = 1;
tau2 = 5;
T = 1*10^3; 

WEE = zeros(N,N);
for i = 1:Q 
    for j = 1:Q
con1 = [mod(i,Q)+1,j];
con2 = [mod(i-2,Q)+1,j];
con3 = [i,mod(j-2,Q)+1];
con4 = [i,mod(j,Q)+1];     
ind = sub2ind([Q,Q],i,j);
indc1 = sub2ind([Q,Q],con1(1),con1(2)); 
indc2 = sub2ind([Q,Q],con2(1),con2(2));
indc3 = sub2ind([Q,Q],con3(1),con3(2)); 
indc4 = sub2ind([Q,Q],con4(1),con4(2));
WEE(ind,indc1)=WE/4;
WEE(ind,indc2)=WE/4;
WEE(ind,indc3)=WE/4;
WEE(ind,indc4)=WE/4;
    end
end



%% integrate and save results
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
save(sprintf('Lattice_sim_WE=%3.3f.mat',WE));


%% plot 
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