title('加入周期信號較小，不足以翻轉')
xlabel('x');
ylabel('U(x)');
v=moviein(20);
subplot(2,1,2)%周期信號較大時
A=3;
[t,yy]=ode23('xinhao1',[0:0.1:10],[1,b-c+A],[],b,z,c,A,w);
n=length(t);
fori=1:fix(n/2)
y=-b.*x.^2/2+c.*x.^4/4-A.*x.*sin(w*t(i));
plot(x,y);
holdon
uu=yy(i,2)*yy(i,2)/2-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
uu=-b.*yy(i,1).^2/2+c.*yy(i,1).^4/4-A.*yy(i,1).*sin(w*t(i));
plot(yy(i,1),uu,'.r','markersize',10)
holdoff
v(:,i)=getframe;




%oula.m子程序:歐拉法解方程
function[x,y,r]=oula(fun,x0,xf,y0,h,r)
n=fix((xf-x0)/h);
x(1)=x0;
y(1)=y0;
x(n)=0;
y(n)=0;
fori=1:(n-1)
x(i+1)=x0+i*h;
y1=y(i)+h*feval(fun,x(i),y(i),r);%調用方程
y2=y(i)+h*feval(fun,x(i+1),y1,r);
y(i+1)=(y1+y2)/2;
end


%f.m方程文件
%參數b,c,A,w為1,故省略。
functionf=f(t,x,r)
f=x-x.^3+cos(t)-r;


%odefilexinhao1.mode文件
functionydot=xinhao1(t,x,flag,b,z,c,A,w)
ydot=[b*x(1)-c*x(1)^3+A*cos(w*t)+z;...
b*x(2)-3*c*x(1)*x(1)*x(2)-A*w*sin(w*t)];