clear

dt=0.005;
tmax = 1.;

Ncells = 50;
dtheta = pi/Ncells;
r=10*ones(1,Ncells);
h=zeros(1,Ncells);
Itot=zeros(1,Ncells);

A = 100.;
epsilon = 0.1;
J0 = -5;
J2 = 5.;
tau = 0.05;

cuecell = 25;

Ith = 10;
Iwidth = 10;
rmax = 100;

for c = 0.1:0.2:0.5;

    for cell = 1:Ncells
        h(cell) = A*c*(1-epsilon + epsilon*cos(2*pi*(cell-cuecell)/Ncells));
    end

    for t = 0.0:dt:tmax
        i = round(t/dt);
        for cell = 1:Ncells
            Itot(cell) = h(cell);
            for i = 1:Ncells
                Itot(cell) = Itot(cell) + r(i)*(J0+J2*cos(2*pi*(i-cell)/Ncells))*dtheta;
            end
        end

        rinf = rmax./(1+exp(-(Itot-Ith)/Iwidth)).^4;
        
        r = rinf + (r-rinf)*exp(-dt/tau);

    end

    hold on
    plot(r)


end

figure
Iin = -100:0.1:100;
rout = rmax./(1+exp(-(Iin-Ith)/Iwidth)).^2;

plot(Iin,rout)