syms u m M L d1 d2 theta dot_theta x dot_x g

f_1 = dot_x ;
f_2 = (u  + m*L*dot_theta^2 * sin(theta) - m*g*cos(theta) *  sin(theta) ) / ( M+m-m* cos(theta)^2 )  - d1*dot_x;
f_3 = dot_theta;
f_4 = (g*sin( theta ) -  cos( theta )*f_2 ) / L - d2*dot_theta;

%% A
var = [x, dot_x, theta, dot_theta];
fun = [f_1, f_2, f_3, f_4];

for i =1 : length(fun)
    for j =1 : length(var)
        disp([i,j])
        diff(fun(i), var(j))
    end
end

%% B
for i =1 : length(fun)
        diff(fun(i), u)
end