function dydt = neuron_Vxx(t,y)
    C = 1;
    Gk = 36;
    GNa = 120;
    Gl = 0.3;
    Ek = -72;
    ENa = 55;
    El = -49.4;
    
    d=0.5;
    N = size(y,1)/4;
    dydt = zeros(N*4,1);
    
    n = y(1);
    m = y(2);
    h = y(3);
    V = y(4);
    
    dydt(1) = (1-y(1))*alphan(y(4))-y(1)*betan(y(4));
    dydt(2) = (1-y(2))*alpham(y(4))-y(2)*betam(y(4));
    dydt(3) = (1-y(3))*alphah(y(4))-y(3)*betah(y(4));
    dydt(4) = -(1/C)*(Gk*y(1)^(4)*(y(4)-Ek)+GNa*y(2)^(3)*y(3)*(y(4)-ENa)+Gl*(y(4)-El))+(y(4*N)-2*V + y(8))/d^2;
    
    for i = 2:N-1
        n = y((i-1)*4+1);
        m = y((i-1)*4+2);
        h = y((i-1)*4+3);
        V = y((i-1)*4+4);
        
        dydt((i-1)*4+1) = ((1-n)*alphan(V))-(n*betan(V));
        dydt((i-1)*4+2) = ((1-m)*alpham(V))-(m*betam(V));
        dydt((i-1)*4+3) = ((1-h)*alphah(V))-(h*betah(V));
        dydt((i-1)*4+4) = -(1/C)*(Gk*n^(4)*(V-Ek)+GNa*m^(3)*h*(V-ENa)+Gl*(V-El))+(y(i*4+4)-2*V + y((i-2)*4+4))/d^2;
    
    end
    
    n = y((N-1)*4+1);
    m = y((N-1)*4+2);
    h = y((N-1)*4+3);
    V = y((N-1)*4+4);
        
    dydt(4*(N-1)+1) = ((1-n)*alphan(V))-(n*betan(V));
    dydt(4*(N-1)+2) = ((1-m)*alpham(V))-(m*betam(V));
    dydt(4*(N-1)+3) = ((1-h)*alphah(V))-(h*betah(V));
    dydt(4*(N-1)+4) = -(1/C)*(Gk*n^(4)*(V-Ek)+GNa*m^(3)*h*(V-ENa)+Gl*(V-El))+(-1*V + y(4*(N-2)+4))/d^2;
    
end