function u = time_varying_LQR(in)

    m = 0.1;
    L = 1;
    g = 9.81;
    M = 1;
    d1 = 1;
    d2 = 0.5;

    % obtain matrices A,B by linearization
    dot_x = in(2);
    theta = in(3);
    dot_theta = in(4);
    x = in(1:4);
    
    Q = diag(in(5:8));
    R = in(9);
    u_ctrl = in(10);
    Ts = in(11);

    f_23 = (g*m*sin(theta)^2 - g*m*cos(theta)^2 + L*dot_theta^2*m*cos(theta))/(- m*cos(theta)^2 + M + m) - (2*m*cos(theta)*sin(theta)*(L*m*sin(theta)*dot_theta^2 + u_ctrl - g*m*cos(theta)*sin(theta)))/(- m*cos(theta)^2 + M + m)^2;
    f_24 = (2*L*dot_theta*m*sin(theta))/(- m*cos(theta)^2 + M + m);
    f_42 =  (d1*cos(theta))/L;
    f_43 =  -(cos(theta)*((g*m*sin(theta)^2 - g*m*cos(theta)^2 + L*dot_theta^2*m*cos(theta))/(- m*cos(theta)^2 + M + m) - (2*m*cos(theta)*sin(theta)*(L*m*sin(theta)*dot_theta^2 + u_ctrl - g*m*cos(theta)*sin(theta)))/(- m*cos(theta)^2 + M + m)^2) + sin(theta)*(d1*dot_x - (L*m*sin(theta)*dot_theta^2 + u_ctrl - g*m*cos(theta)*sin(theta))/(- m*cos(theta)^2 + M + m)) - g*cos(theta))/L;
    f_44 =  - d2 - (2*dot_theta*m*cos(theta)*sin(theta))/(- m*cos(theta)^2 + M + m);

    A = [0,1,0,0 ; 0,-d1,f_23, f_24; 0,0,0,1; 0, f_42, f_43, f_44 ] ;

    B = [0;
        1/(- m*cos(theta)^2 + M + m);
        0 ;
        -cos(theta)/(L*(- m*cos(theta)^2 + M + m))  ]; 

    Ad = expm(A*Ts);
    Bd = pinv(A)*(Ad-eye(4))*B;
    
    try
        [X,L,K_lqr_d] = dare(Ad,Bd,Q, R);
        u = -K_lqr_d*x;
    catch
        u = -1000;
    end
    

end
%time_varying_LQR