1、首先，获取系统模型方程。
2、其次，调用dsolve 函数求得方程的解：
>> s = dsolve('2*D2x+Dx+8*x = 0','Dp = -p - 17*x','x(0)=4',' Dx(0)=0','p(0)=0')
3、接着，访问方程的解：
>> s.x
ans =
exp(-1/4*t)*(4/21*sin(3/4*7^(1/2)*t)*7^(1/2)+4*cos(3/4*7^(1/2)
*t))
>> s.p
ans =
-68/9*cos(3/4*7^(1/2)*t)*exp(-1/4*t)-748/63*7^(1/2)*exp(-1/4*t
)*sin(3/4*7^(1/2)*t)+68/9*exp(-t)
4、绘制位置随时间变化的函数图象：
>> ezplot(s.x,[0 10])
>> title('质量的位置')
5、绘制在相同时间内的动量：
>> ezplot(s.p,[0 10])
>> title('动量')
6、调用plot 产生相图：
>> plot(xval,pval),xlabel('x'),ylabel('p'),title('质量弹簧的相图')
7、代码总结：
>> s = dsolve('2*D2x+Dx+8*x = 0','Dp = -p - 17*x','x(0)=4',' Dx(0)=0','p(0)=0')；
>> s.x
ans =
exp(-1/4*t)*(4/21*sin(3/4*7^(1/2)*t)*7^(1/2)+4*cos(3/4*7^(1/2)
*t))
>> s.p
ans =
-68/9*cos(3/4*7^(1/2)*t)*exp(-1/4*t)-748/63*7^(1/2)*exp(-1/4*t
)*sin(3/4*7^(1/2)*t)+68/9*exp(-t)
>> ezplot(s.x,[0 10])
>> title('质量的位置')
>> ezplot(s.p,[0 10])
>> title('动量')
>> ezplot(s.x,s.p,[-5 5])
>> axis([-8 8 -25 20])
>> title('参量图')
>> tvalues = (0:0.1:10);
>> xval = subs(s.x,'t',tvalues);
>> pval = subs(s.p,'t',tvalues);
>> plot(xval,pval),xlabel('x'),ylabel('p'),title('质量弹簧的相图')
