function d_theta = ode_simulation_kuka(t,theta)

    %disp(t)

    global l_1 l_2 k1 k2 p_i p_1_0 p_2_0 gamma_1 gamma_2 K step_size_angle step_size_velocity
    global q_prev_theta_1 q_prev_theta_2 q_prev_theta_3 q_prev_theta_4
    global k_1 k_2 % Joints Friction Coeefficients
    global I_1 I_2 m_1 m_2 l1 l2 g

    % Matrices Computations
    M_11 = I_1 + I_2 + m_1*((l1^2)/4) + m_2*((l1^2) + ((l2^2)/4) + l1*l2*cos(theta(3)));
    M_12 = I_2 + m_2*(((l2^2)/4) + (l1*l2*cos(theta(3)))/2);
    M_21 = M_12;
    M_22 = I_2 + m_2*((l2^2)/4);
    M = [M_11 M_12 ; M_21 M_22];
    c = (m_2*l1*l2*sin(theta(3)))/2;
    C = [-c*theta(4) -c*(theta(2,:)+theta(4)) ; ...
        c*theta(2) 0];
    G = [((m_1*g*l1*cos(theta(1)))/2)+(m_2*g*(l1*cos(theta(1))+(l2*cos(theta(1)+theta(3)))/2)) ; ...
        (m_2*g*l2*cos(theta(1)+theta(3)))/2];

    % Quantize Measurements
    if t ~= 0
        q_theta_1 = uniform_hysteretic_quantizer(theta(1), q_prev_theta_1, step_size_angle/2);
        q_theta_2 = uniform_hysteretic_quantizer(theta(2), q_prev_theta_2, step_size_velocity/2);
        q_theta_3 = uniform_hysteretic_quantizer(theta(3), q_prev_theta_3, step_size_angle/2);
        q_theta_4 = uniform_hysteretic_quantizer(theta(4), q_prev_theta_4, step_size_velocity/2);
        q_theta = [q_theta_1 q_theta_2 q_theta_3 q_theta_4];
    else
        q_theta_1 = q_prev_theta_1;
        q_theta_2 = q_prev_theta_2;
        q_theta_3 = q_prev_theta_3;
        q_theta_4 = q_prev_theta_4;
        q_theta = [q_theta_1 q_theta_2 q_theta_3 q_theta_4];
    end
    q_prev_theta_1 = q_theta_1;
    q_prev_theta_2 = q_theta_2;
    q_prev_theta_3 = q_theta_3;
    q_prev_theta_4 = q_theta_4;

    % Desired Angle Trajectories
    if t <= 6

        theta_d_1 = 0.01454441*t^3 - 0.0036361*t^4 + 0.00024241*t^5;
        theta_dot_d_1 = 3*0.01454441*t^2 - 4*0.0036361*t^3 + 5*0.00024241*t^4;
        %theta_d_1 = 0.5*sin(0.3*pi*t);
        %theta_dot_d_1 = 0.5*0.3*pi*cos(0.3*pi*t);

        theta_d_2 = 0.01454441*t^3 - 0.0036361*t^4 + 0.00024241*t^5;
        theta_dot_d_2 = 3*0.01454441*t^2 - 4*0.0036361*t^3 + 5*0.00024241*t^4;
        %theta_d_2 = 0.5*cos(0.3*pi*t);
        %theta_dot_d_2 = -0.5*0.3*pi*sin(0.3*pi*t);

    else
        theta_d_1 = pi/10;
        theta_dot_d_1 = 0;

        theta_d_2 = pi/10;
        theta_dot_d_2 = 0;

    end

    % Disturbance Signals
    d_1 = 0.25*sin(t);
    d_2 = 0.75*cos(t);

    % Performance Functions --> Define the envelope in which the output
    % (tracking) error resides
    p_1 = p_1_0*exp(-l_1*t)+k1*p_i;
    p_2 = p_2_0*exp(-l_2*t)+k2*p_i;

    % Construct Essential Signal for the final Control Signal
    delta_1 = k1*step_size_velocity+k1*step_size_angle;
    delta_2 = k2*step_size_velocity+k2*step_size_angle;

    e_q_1 = q_theta(1)-theta_d_1;
    dot_e_q_1 = q_theta(2)-theta_dot_d_1;
    s_1_q = (dot_e_q_1+k1*e_q_1)/(p_1-delta_1);
    T_dot_s_q_1 = -(2/((s_1_q^2)-1));

    e_q_2 = q_theta(3)-theta_d_2;
    dot_e_q_2 = q_theta(4)-theta_dot_d_2;
    s_2_q = (dot_e_q_2+k2*e_q_2)/(p_2-delta_2);
    T_dot_s_q_2 = -(2/((s_2_q^2)-1));

    u_1 = (K*gamma_1*(p_1-delta_1)*s_1_q)/(T_dot_s_q_1*(gamma_1*(s_1_q^2)+gamma_2*(s_2_q^2)-1));
    u_2 = (K*gamma_2*(p_2-delta_2)*s_2_q)/(T_dot_s_q_2*(gamma_1*(s_1_q^2)+gamma_2*(s_2_q^2)-1));
    u = [u_1 ; u_2];

    inv_M = inv(M);

    d_theta =[ theta(2) ; ...
        -inv_M(1,:)*C*[theta(2) ; theta(4)]-inv_M(1,:)*G+inv_M(1,:)*u+d_1 ;  ...
        theta(4) ; ...
        -inv_M(2,:)*C*[theta(2) ; theta(4)]-inv_M(2,:)*G+inv_M(2,:)*u+d_2 ...
        ];

end