function xdot=molla(t,x,flag,u)

global beta k1 k2 M1 M2

xdot=zeros(4,1);

xdot(1)=x(2);
xdot(2)=-(k1+k2)/M1*x(1)-beta/M1*x(2)+k2/M1*x(3)+u(2);
xdot(3)=x(4);
xdot(4)=k2/M2*x(1)-k2/M2*x(3)+u(4);