这个程序怎么改啊？？是关于陀螺反馈线性化的

A=[0 1 0 0;a21 a22 a23 a24;0 0 0 1;a41 a42 a43 a44];B=[0 0;1 0;0 0;0 1];C=[1 0 0 0;0 0 1 0];ommiga_1=6635.6;ommiga_2=6632.3;tau_1=0.8;tau_2=0.9;theta_ommiga= sind(10);theta_tau=sind(5);k=0.4;n=3;Ommiga=0;a21=-1\ommiga^2+ommigaDeltaommiga*cos(2\theta_ommiga);a22=-2/tau-Delta(1/tau)*cos(2\theta_tau);a23=ommigaDeltaommiga*sin(2\theta_ommiga);a24=2*k*Ommiga-Delta(1/tau)*sin(2\theta_tau);a41=ommigaDeltaommiga*sin(2\theta_ommiga);a42=2*k*Ommiga-Delta(1/tau)*sin(2\theta_tau);a43=-1\ommiga^2-ommigaDeltaommiga*cos(2\theta_ommiga);a44=-2/tau+Delta(1/tau)*cos(2\theta_\tau);ommiga^2==((ommiga_1)^2+(ommiga_2)^2)/2;ommigaDeltaommiga==((ommiga_1)^2-(ommiga_2)^2)/2;1/tau==1/2*(1/tau_1+1/tau_2);Delta(1/tau)==1/2*(1/tau_1-1/tau_2);sys=ss(A,B,C);t=[0:0.01:20];u=[f_x;f_y];f_x=x"-2*n*k*Ommiga*y'+(2/tau+Delta(1/tau)*cos(2*n*theta_tau))*x'+(Delta(1/tau)*sin(2*n*theta_tau))*y'+(ommiga^2-ommigaDeltaommiga*cos(2*n*theta_ommiga)*x-ommigaDeltaommiga*sin(2*n*theta_ommiga)*y;f_y=y"+2*n*k*Ommiga*x'+(2/tau-Delta(1/tau)*cos(2*n*theta_tau))*y'+(Delta(1/tau)*sin(2*n*theta_tau))*x'+(ommiga^2+ommigaDeltaommiga*cos(2*n*theta_ommiga)*y-ommigaDeltaommiga*sin(2*n*\theta_\ommiga)*x;x0=[0;0];[y,T,x]=lsim(sys,u,t,x0);subplot(2,1,1),plot(T,x(:,1));xlabel('Time(sec)'),ylabel('output');