如何利用一下程序作出poincare图或分叉图

是滑动轴承在某一参数下的动力学特性：

function df=rhs(t,f,data)

%方程中的系统参数

m=data(1);g=data(2);c=data(3);n=data(4);b=data(5);mf=data(6);M=data(7);D0=data(8);De=data(9);

tau0=data(10);Ke=data(11);K0=data(12);omega=data(13);e=data(14);

%化简后的动力学方程

df(1)=f(2);

df(2)=-(De+D0*(1-f(1)^2-f(3)^2)^(-n))*f(2)/(M*omega)-2*tau0*mf*f(4)*(1-sqrt(f(1)^2+f(3)^2))^b/M-(Ke+K0*(1-f(1)^2-f(3)^2)^(-n)-mf*omega^2*tau0^2*(1-sqrt(f(1)^2+f(3)^2))^(2*b))*f(1)/(M*omega^2)-(tau0*D0*(1-sqrt(f(1)^2+f(3)^2))^b*(1-f(1)^2-f(3)^2)^(-n))*f(3)/(M*omega)+m*e*cos(t)/(M*c);

df(3)=f(4);

df(4)=-(De+D0*(1-f(1)^2-f(3)^2)^(-n))*f(4)/(M*omega)+2*tau0*mf*f(2)*(1-sqrt(f(1)^2+f(3)^2))^b/M-(Ke+K0*(1-f(1)^2-f(3)^2)^(-n)-mf*omega^2*tau0^2*(1-sqrt(f(1)^2+f(3)^2))^(2*b))*f(3)/(M*omega^2)+(tau0*D0*(1-sqrt(f(1)^2+f(3)^2))^b*(1-f(1)^2-f(3)^2)^(-n))*f(1)/(M*omega)+m*e*sin(t)/(M*c)-m*g/(M*c*omega^2);

df=[df(1);df(2);df(3);df(4)];

运行以下程序

clc;

clear;

format long

m=50;       Ke=7.2762e6;    De=2e3;

n0=0.079;    m0=-0.25;    g=9.8;

e=0.0002;    R=0.067;      l=0.102;

n=2.5;       b=0.45;     tau0=0.4;   c=0.0003;      nu=50;

z=0.1;       p=2e5;        delta=0.0003;    T=l/nu;     upsilon=1.47e-5;

fid=fopen('poincare.txt','w');

fidd=fopen('time.txt','w');

fiddd=fopen('center.txt','w');

% 转速

omega=1300;

% 初始位置

x0=0;

y0=-m*g/Ke;

X0=x0/c;Y0=y0/c;

% 无量纲转速

nq=length(omega);

for i=1:nq

s(i)=omega(i)/sqrt(Ke/m);

% Child短密封轴承动力系数

Rv(i)=omega(i)*R*delta/upsilon;

Ra=2*nu*delta/upsilon;

lambda(i)=n0*Ra^m0*(1+(Rv(i)/Ra)^2)^((1+m0)/2);

sigma(i)=lambda(i)*l/delta;

B(i)=2-((Rv(i)/Ra)^2-m0)/((Rv(i)/Ra)^2+1);

E(i)=(1+z)/2/(1+z+2*sigma(i));

u0(i)=2*sigma(i)^2*E(i)*(1-m0)/(1+z+2*sigma(i));

u1(i)=2*sigma(i)^2*(E(i)/sigma(i)+B(i)*(1/6+E(i))/2)/(1+z+2*sigma(i));

u2(i)=sigma(i)*(1/6+E(i))/(1+z+2*sigma(i));

u3(i)=pi*R*p/lambda(i);

% 密封力Muszynska模型系数

% K=K0*(1-e^2)^(-n)

% D=D0*(1-e^2)^(-n)

% tau=tau0*(1-e)^b

% e=sqrt(x^2+y^2)/c

K0(i)=u0(i)*u3(i);

D0(i)=u1(i)*u3(i)*T;

mf(i)=u2(i)*u3(i)*(T^2);

M(i)=mf(i)+m;

% 系统结构动力学方程中的参数

data=[m,g,c,n,b,mf(i),M(i),D0(i),De,tau0,Ke,K0(i),omega(i),e];

%初始条件

f0=[X0;0;Y0;0];

t0=0;h=0.001;tf=600;

tn=[t0:h:tf];

[t,f]=ode45(@rhs,tn,f0,[],data);

%求出的是无量纲化的位移、速度

%现转化为实际的位移、速度

x=f(:,1);y=f(:,3);

dx=f(:,2);dy=f(:,4);

nn=(tf-t0)/h+1;nnn=300000;%nnn为除去瞬态响应后的第一点数，nn为采集的总点数

nnnn=400000;

for i=nnn:nn

fprintf(fidd,'%6.4f  %12.8f\n',t(i),x(i));hold on;

end

fclose(fidd);

for j=nnn:nn

fprintf(fiddd,'%12.8f %12.8f\n',x(j),y(j));hold on;

end

fclose(fiddd);

% % 求庞加莱映射点

[A,jj]=max(x(nnn:nnnn));

jj=jj+nnn-1;

TT=2*pi;

for k=1:nn

tt(k)=(jj-1)*h+k*TT;

if tt(k)>tf

break;

end

for kk=nnn:nn

if tt(k)>(kk-1)*h&tt(k)

l1=(tt(k)-(kk+1)*h)*(tt(k)-kk*h)/(2*h^2);

l2=(tt(k)-(kk-1)*h)*(tt(k)-(kk+1)*h)/(-h^2);

l3=(tt(k)-(kk-1)*h)*(tt(k)-kk*h)/(2*h^2);

xx(k)=l1*x(kk-1)+l2*x(kk)+l3*x(kk+1);

dx(k)=l1*dx(kk-1)+l2*dx(kk)+l3*dx(kk+1);

fprintf(fid,'%12.8f %12.8f\n',xx(k),dx(k));hold on;

end

end

end

end

fclose(fid);