function [theta_store,centroidStore] = walking_diff_eq_function(numStrides,h,deltaVec,p,phiVec,c_2_vec,omega,visualization_flag)

% NOTE: switch between 5 and 6 legs just set c_2_vec(4) = 0 (turn off front right
% leg force or set delta_vec(4) = 0, never put it down). The former is preferable.

    M = 1; % mass
    B1 = 1.5; % linear damping
    I = .01; % moment of inertia
    B2 = 1.5; % angular damping

    
    L = 1;
    
    % 
    % %%%% Legs
    tau_L = 0.01; % leg relaxation timescale
    c_1_vec = [0 0 0 0 0 0]; % leg force perpendicular to body axis
    l1_star = 0.59; % relaxed length leg i
    theta1_star = 1.19; % relaxed angle leg i
    l2_star = 0.66;
    theta2_star = 0.33;
    l3_star = 0.42;
    theta3_star = -.69;
    l4_star = 0.59;
    theta4_star = 1.19;
    l5_star = 0.66;
    theta5_star = 0.33;
    l6_star = 0.42;
    theta6_star = -0.69;
    p1_x = -0.05; % leg attachment to body
    p1_y = 0.2;
    p2_x = -0.05;
    p2_y = 0;
    p3_x = -0.05; 
    p3_y = -0.11;
    p4_x = 0.05; 
    p4_y = 0.2;
    p5_x = 0.05; 
    p5_y = 0; 
    p6_x = 0.05; 
    p6_y = -0.11;


    n_hat = 0.9; % neural threshold for firing

    %%%% Neurons
    tau_n = .01; % neuron relaxation timescale (fast (small) compared to min time neuron on by CPG)

    v_I = 1;
    alpha = 1; % strength of excitation pulse, has to be greater than n_hat



    X_init = 0;
    X_dot_init = 0;
    Y_init = 0;
    Y_dot_init = 0;
    Theta_init = 0;
    Theta_dot_init = 0;

    bodyInitConditions = [X_init,X_dot_init,Y_init,Y_dot_init,Theta_init,Theta_dot_init];

    %%%% Legs
    x1_init = -l1_star*cos(theta1_star)+p1_x+X_init;
    y1_init = l1_star*sin(theta1_star)+p1_y+Y_init;
    x2_init = -l2_star*cos(theta2_star)+p2_x+X_init;
    y2_init = l2_star*sin(theta2_star)+p2_y+Y_init;
    x3_init = -l3_star*cos(theta3_star)+p3_x+X_init;
    y3_init = l3_star*sin(theta3_star)+p3_y+Y_init;
    x4_init = l4_star*cos(theta4_star)+p4_x+X_init;
    y4_init = l4_star*sin(theta4_star)+p4_y+Y_init;
    x5_init = l5_star*cos(theta5_star)+p5_x+X_init;
    y5_init = l5_star*sin(theta5_star)+p5_y+Y_init;
    x6_init = l6_star*cos(theta6_star)+p6_x+X_init;
    y6_init = l6_star*sin(theta6_star)+p6_y+Y_init;


    leg_pos_matrix_init = [x1_init,x2_init,x3_init,x4_init,x5_init,x6_init;
        y1_init,y2_init,y3_init,y4_init,y5_init,y6_init];

    legInitConditions = reshape(leg_pos_matrix_init,1,12);
 
    %%%% Neurons
    n1_init = 0;
    n2_init = 0;
    n3_init = 0;
    n4_init = 0;
    n5_init = 0;
    n6_init = 0;
    neuroInitConditions = [n1_init,n2_init,n3_init,n4_init,n5_init,n6_init]; % start with 


    % Final initial conditions vector
    initialConditionsVector = [bodyInitConditions,legInitConditions,neuroInitConditions];

    t_init = 0.00001; % small adjustment from 0, otherwise H(.) evaluate to 1/2, creating some weird effects
    t_final = numStrides/omega;
    tspan = t_init:h:t_final;
    T = length(tspan);
    theta_store = zeros(1,T);
    theta_store(1) = initialConditionsVector(5);
    centroidStore = zeros(2,T);
    centroidStore(:,1) = [initialConditionsVector(1),initialConditionsVector(3)];

    num_equations = length(initialConditionsVector); 
    sol = zeros(num_equations,T);
    
    sol(:,1) = initialConditionsVector;
    
    for i = 2:length(tspan)

        t_i = tspan(i-1);
        y_old = sol(:,i-1);

        X = y_old(1);
        X_dot = y_old(2);
        Y = y_old(3);
        Y_dot = y_old(4);
        Theta = y_old(5);
        Theta_dot = y_old(6);

        %%%% Legs
        x1 = y_old(7);
        y1 = y_old(8);
        x2 = y_old(9);
        y2 = y_old(10);
        x3 = y_old(11);
        y3 = y_old(12);
        x4 = y_old(13);
        y4 = y_old(14);
        x5 = y_old(15);
        y5 = y_old(16);
        x6 = y_old(17);
        y6 = y_old(18);

        %%%% Neurons
        n1 = y_old(19);
        n2 = y_old(20);
        n3 = y_old(21);
        n4 = y_old(22);
        n5 = y_old(23);
        n6 = y_old(24);

        R_Theta = [cos(Theta),sin(Theta);-sin(Theta),cos(Theta)]; % rotation to lab frame


        % neuron sigmoid calculation (on/off)

        neuron_sigmoid_vec = [heaviside(n1-n_hat),...
                heaviside(n2-n_hat),...
                heaviside(n3-n_hat),...
                heaviside(n4-n_hat),...
                heaviside(n5-n_hat),...
                heaviside(n6-n_hat)];

        % force calculation
        force_legs = R_Theta*[c_1_vec;c_2_vec].*repmat(neuron_sigmoid_vec,2,1);
            
        leg_pos_vec = [x1,x2,x3,x4,x5,x6;y1,y2,y3,y4,y5,y6];
        r_vec = leg_pos_vec-repmat([X;Y],1,6); % distance from leg endpoint to COM

        F_vec_temp = 1/M*(sum(force_legs,2) - B1*[X_dot;Y_dot]);
        leg_torque_vec = r_vec(1,:).*force_legs(2,:)-r_vec(2,:).*force_legs(1,:);
        m_temp = 1/I*(sum(leg_torque_vec) - B2*Theta_dot);

        p_temp = [p1_x p2_x p3_x p4_x p5_x p6_x;p1_y p2_y p3_y p4_y p5_y p6_y];
        relaxed_temp = R_Theta*(leg_pos_matrix_init) + repmat([X;Y],1,6);

        leg_dynamics_temp = 1/tau_L*(relaxed_temp-leg_pos_vec).*repmat(ones(1,6)-neuron_sigmoid_vec,2,1); % firing = exert force (stay in place), not firing = relax

        neuron_vec = [n1 n2 n3 n4 n5 n6];

        leg_vec_from_pivot = leg_pos_vec - (R_Theta*p_temp + repmat([X;Y],1,6));
        leg_length_vec = diag(sqrt(leg_vec_from_pivot'*leg_vec_from_pivot))';
        leg_sigmoid_vec_temp = heaviside(leg_length_vec-[(1+p)*l1_star,(1+p)*l2_star,(1+p)*l3_star,(1+p)*l4_star,(1+p)*l5_star,(1+p)*l6_star]) +...
            heaviside([(1-p)*l1_star,(1-p)*l2_star,(1-p)*l3_star,(1-p)*l4_star,(1-p)*l5_star,(1-p)*l6_star]-leg_length_vec); % (1+p)*li_star, (1-p)*li_star bounds
                % for extension and compression. So, if a leg is too stretched or too compressed, this evaluates to one for that leg - subtract from firing.

        neuron_equation_vec = 1/tau_n*(-neuron_vec - v_I*leg_sigmoid_vec_temp + alpha*(mod(omega*t_i*ones(1,6)-phiVec,1)>(1-deltaVec))); % default state: off, except modulated by alpha*[T/F criterion]
          
        %%%%%% BASIC VISUALIZATION
    if visualization_flag == 1
        if mod(i,1/10*1/h) == 0 % note: higher h will make things seem to move slower due to more calculating time between frames
            subplot(2,1,1)
            plot(X,Y,'ko','MarkerSize',6)
            axis([-5,5,-5,10])
            hold on

            plot(leg_pos_vec(1,:),leg_pos_vec(2,:),'o','MarkerSize',6)

            p1 = R_Theta*[p1_x;p1_y]+[X;Y];
            p2 = R_Theta*[p2_x;p2_y]+[X;Y];
            p3 = R_Theta*[p3_x;p3_y]+[X;Y];
            p4 = R_Theta*[p4_x;p4_y]+[X;Y];
            p5 = R_Theta*[p5_x;p5_y]+[X;Y];    
            p6 = R_Theta*[p6_x;p6_y]+[X;Y];

            plot(p1(1),p1(2),'ko','MarkerSize',5)
            plot(p2(1),p2(2),'ko','MarkerSize',5)
            plot(p3(1),p3(2),'ko','MarkerSize',5)
            plot(p4(1),p4(2),'ko','MarkerSize',5)
            plot(p5(1),p5(2),'ko','MarkerSize',5)
            plot(p6(1),p6(2),'ko','MarkerSize',5)


            hold off

            subplot(2,1,2)
            hold on
            plot(t_i,neuron_vec(1),'r+','MarkerSize',5)
            plot(t_i,neuron_vec(2),'g+','MarkerSize',5)
            plot(t_i,neuron_vec(3),'b+','MarkerSize',5)
            plot(t_i,neuron_vec(4),'ro','MarkerSize',5)
            plot(t_i,neuron_vec(5),'go','MarkerSize',5)
            plot(t_i,neuron_vec(6),'bo','MarkerSize',5)
            axis([0,5,0,1])

            hold off
            pause(.1)
        end
    end


        %%%%%% END VISUALIZATION    
        

        % EULER'S METHOD
        dy(1) = X_dot;
        dy(2) = F_vec_temp(1);
        dy(3) = Y_dot;
        dy(4) = F_vec_temp(2);
        dy(5) = Theta_dot;
        dy(6) = m_temp;

        %%%% Legs

        dy(7) = leg_dynamics_temp(1,1);
        dy(8) = leg_dynamics_temp(2,1);
        dy(9) = leg_dynamics_temp(1,2);
        dy(10) = leg_dynamics_temp(2,2);
        dy(11) = leg_dynamics_temp(1,3);
        dy(12) = leg_dynamics_temp(2,3);
        dy(13) = leg_dynamics_temp(1,4);
        dy(14) = leg_dynamics_temp(2,4);
        dy(15) = leg_dynamics_temp(1,5);
        dy(16) = leg_dynamics_temp(2,5);
        dy(17) = leg_dynamics_temp(1,6);
        dy(18) = leg_dynamics_temp(2,6);

        %%%% Neurons
        dy(19:24) = neuron_equation_vec;

        sol(:,i) = y_old + h*dy';

        current_solution = sol(:,i);
        theta_store(i) = current_solution(5);
        centroidStore(:,i) = [current_solution(1),current_solution(3)];
    end

end