%----------------------------------------------------------
% 專門用于求解病態系統的算法
%----------------------------------------------------------
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[t,Y] = ode15s(@rigid,[0:0.01:1.5],[1 0 -1],options);

%----------------------------------------------------------
%解析解
%----------------------------------------------------------
x1=1/2*exp(-40*t).*cos(40*t)+1/2*exp(-2*t)+1/2*exp(-40*t).*sin(40*t) ;
x2=-1/2*exp(-40*t).*cos(40*t)+1/2*exp(-2*t)-1/2*exp(-40*t).*sin(40*t);
x3=exp(-40*t).*sin(40*t)-exp(-40*t).*cos(40*t);

%----------------------------------------------------------
%通過輸出圖像進行比較 
%----------------------------------------------------------
subplot(3,1,1)
plot(t,x1-Y(:,1),'.r'),xlabel('t'),ylabel('y1')           

subplot(3,1,2)
plot(t,x2-Y(:,2),'.r'),xlabel('t'),ylabel('y2')

subplot(3,1,3)
plot(t,x3-Y(:,3),'.r'),xlabel('t'),ylabel('y3')