%図7.1, 図7.2
clear all; close all;

%システムパラメータを与える
A = [0 1;-6 -5]; %行列A
b = [0;1]; %ベクトルb
c = [1 0;0 1]; %x_1とx_2をプロットするためにcを単位行列にする
d = [0;0]; %cに合わせてdは零ベクトルとする

%システムの状態空間表現を与える
sysP = ss(A,b,c,d); %状態空間表現

%システム行列Aの固有値を求める
eig(A) %システム行列Aの固有値

%時間変数の定義
t = 0:0.01:10; %0から10まで0.01刻み

%ステップ信号の生成(初期値が異なるステップ応答を求めるため)
t1 = size(t); tt = t1(:,2); %時間変数のサイズを抽出
u = ones(1,tt); %時間変数のサイズに応じたステップ信号の生成

%システムの初期値を与える
x0 = [1;0]; %x(0) = [1;0]'の場合
x00 = [0;0]; %x(0) = [0;0]'の場合

%ステップ応答の計算
y = lsim(sysP,u,t,x0); %x(0) = [1;0]'の場合のステップ応答
y_0 = lsim(sysP,u,t,x00); %x(0) = [0;0]'の場合のステップ応答

figure(1) %図のウンドウを開く
plot(t,y(:,1),'-b',t,y(:,2),'-r'); %x_1とx_2を抽出してプロット
xlim([0 10]); %横軸(時間軸)の範囲の指定
ylim([-1.0 1.0]); %縦軸の範囲の設定
grid; %罫線で表示
xlabel('time t[s]'); %横軸のラベル表示
ylabel('y(t)'); %縦軸のラベル表示
legend('x_1(t)','x_2(t)'); %凡例の表示

figure(2) %図のウィンドウを開く
plot(t,y_0(:,1),'-b',t,y_0(:,2),'-r'); %x_1とx_2を抽出してプロット
xlim([0 10]); %横軸(時間軸)の範囲の指定
ylim([-1.0 1.0]); %縦軸の範囲の設定
grid; %罫線を表示
xlabel('time t[s]'); %横軸のラベル表示
ylabel('x_1(t),x_2(t)'); %縦軸のラベル表示
legend('x_1(t)','x_2(t)'); %凡例の表示
ans =

   -2.0000
   -3.0000