分岔图分析程序。用于分析转子动力学中振动稳定特性

				非线性微分方程组】怎么在Matlab求解后画分岔图
				
				
				这个是一个Matlab求解非线性方程的程序，不明白求解完后怎么画分岔图，还有庞加莱图，求助高手指点，小弟不胜感激
				
				主程序
				
				[Copy to clipboard] [ - ]CODE:
				
				function BallBrg_NonL_Forum
				% 求解外圈固定球轴承的变柔度(VC-Varying Compliance)振动（基于赵凌燕的论文）
				% 程序有一些不合理、甚至错误的地方，可以用更好的代码代替，由于时间关系没有修改，
				% 如有人感兴趣可以把修改的程序发布出来。
				%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
				% 作者：toes
				% 版本：论坛发布版
				% 相关程序：BallBrg_NonL_Sub_Forum
				% 调试环境：Matlab7.0 WinXP SP2
				% 参考文献：
				% 1.赵凌燕.滚动轴承－转子系统的非线性动力学研究.西北工业大学硕士论文.2003.3.
				%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
				clear
				clc
				
				%% 参数设置
				% 用了全局变量来传递一些变量，不推荐，但是懒得改了，好心人优化一下。
				global w d D Nb gama kn M C F
				
				% 为了方便绘制分岔图而设置的参数
				n_One_T = 100;% 每个周期的采样点数
				n_T = 100;% 采样时间占几个周期
				
				% 61903/P5(17*30*7) 球轴承参数
				d=0.0173;% 内滚道直径
				D=0.0265;% 外滚道直径
				Nb=9; % 滚子数
				
				n_n = 0;
				w_limit1=100;% 最低转速（rpm）
				w_limit2=20000;% 最高转速（rpm）
				w_step = 100;% 转速变化步长(rpm)
				
				q_initial(1:4,1) = 1e-11;% 初始值
				gama = 0.00002;% 间隙(m)
				F = 6;% 径向力(N)
				kn = 7.055e9*0.001^1.5;
				% 滚子与滚道之间接触力与变形量的关系(N/mm^1.5)。赵的论文给出。
				M=0.6*[1 0;0 1];% 质量矩阵
				C=200*[1 0; 0 1];% 阻尼矩阵
				
				%% 响应计算循环
				for w_rpm=w_limit1:w_step:w_limit2
				
				n_n = n_n+1 % 计数变量
				disp(w_rpm)
				w = w_rpm*pi/30;% 转化为rad/s单位
				
				wi = w;% 内圈角速度
				wo = 0;% 外圈角速度
				
				w_cage = ( wi*d/2+wo*D/2 )/2/((D+d)/4);% 保持架
				w_vc = w_cage*Nb/2/pi; % 变刚度频率(vc频率)。单位Hz
				T_vc = 1/w_vc;% vc周期
				
				dt=T_vc/n_One_T;% 取点时间步长，dt随转速变化。
				time=n_T*T_vc;% 总的时间
				
				n = round(time/dt);% 离散点数
				t_span(1:n) = linspace(0,time,n);% 时间数组
				
				[t,q]= ode23tb('BallBrg_NonL_Sub_Forum', t_span, q_initial);
				% 至于用什么ode函数求解合适需要比较验证
				
				Q{n_n}=q;
				save Q.mat Q; % 存储数据
				end
				disp('Calculation is done!')
				
				子程序
				
				[Copy to clipboard] [ - ]CODE:
				
				function dq = BallBrg_NonL_Sub_Forum(t,q)
				% BallBrg_NonL调用的微分方程子程序
				% 求解外圈固定球轴承的变柔度(VC)振动
				%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
				% 作者：toes
				% 版本：论坛发布版
				% 相关程序：BallBrg_NonL_Forum
				% 参考文献：
				% 1.赵凌燕.滚动轴承－转子系统的非线性动力学研究.西北工业大学硕士论文.2003.3.
				%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
				global w d D Nb gama kn M C F
				
				wi = w;
				wo = 0;
				w_cage=( wi*d+wo*D )/4/((D+d)/4);% 保持架转速（rad/s）
				
				fq=zeros(2,1);% 轴承力初值
				
				diff_1_3 = q(1,1);% 水平方向位移
				diff_2_4 = q(2,1);% 垂直方向位移
				
				% 求轴承的非线性反力
				for No_ball=1:Nb
				sita(No_ball) = 2*pi/Nb*(No_ball-1) + w_cage*t;% 第No_ball个滚珠的位置角
				Clearance(No_ball,1) = diff_1_3*sin( sita(No_ball) ) ...
				+ diff_2_4*cos( sita(No_ball) ) - gama;% 滚珠与内滚道的间隙变化。
				% 判断哪几个滚动体受到接触力
				if Clearance(No_ball)				Clearance(No_ball) = 0;
				end
				fs = abs( (1000*Clearance(No_ball))^1.5 );
				
				fq(1,1) = fq(1,1)+kn*fs*sin(sita(No_ball));
				fq(2,1) = fq(2,1)+kn*fs*cos(sita(No_ball));
				end
				
				F_m1d1_cos = 0;% 不平衡力在水平方向的投影。本例不考虑。
				F_m1d1_sin = 0;% 不平衡力在垂直方向的投影。本例不考虑。
				
				Fq(1,1)= - fq(1,1) + F_m1d1_cos;% 水平方向外力
				Fq(2,1)= - fq(2,1) + F_m1d1_sin - F;% 垂直方向外力
				
				K = [0 0; 0 0];% 刚性转子，轴段为刚性。
				
				% 动力学微分方程
				dq(3:4,1)=inv(M)*(Fq-K*q(1:2,1)-C*q(3:4,1));% x和y方向加速度
				dq(1:2,1)=q(3:4,1);