function xdot=vdp(t,x)
%求著名的Van Der Pol 方程x"+( x^2 −1)x' + x = 0的数值解并绘制其时间响应
%曲线和状态轨迹图
%1.演化为状态方程
%令x1 = x',x2= x, 把x"+( x^2 −1)x' + x = 0写成状态方程x1'=(1-x2^2)*x1-x2,x2'=x1
xdot=zeros(2,1);%使xdot 成为二元零向量采用列向量以便被matlab 其他指令调用
xdot(1)=(1-x(2)^2)*x(1)-x(2);
xdot(2)=x(1);
将上面语句保存为vdp.m后,在窗口中输入下面语句看结果：
t0=0;tf=20;x0=[0,0.25]';[t,x]=ode45('vdp',[t0,tf],x0);subplot(1,2,1),plot(t,x(:,1),':b',t,x(:,2),'-r'),
legend('速度','位移'),subplot(1,2,2),plot(x(:,1),x(:,2))