function selkov_parameter_dependency()
    % selcovモデル
    % ODEを数値的に解く
    %% パラメータの指定
    s0 = [1, 1]; % ADP(x)の初期値とF6P(y)の初期値
    param = [0.06, 0.6]; % パラメータ a, b
    time = 0.01:0.05:400; %シミュレーションを行う時間
    b = 0.1:0.2:1.1; % パラメータを振る
    col = jet; % 色指定に使用
    figure(1);
    clf;
    figure(2);
    clf;
    
    %% パラメータを振り分けて計算
    for i = 1:numel(b);
        param(2) = b(i); %% 振ったパラメータを代入
        
        %% 数値計算
        [t, time_course] = ode45(@(t, s) ODE(t, s, param), time, s0); % ODEを解く
        
        %% figure 1 を描く
        figure(1);
        hold on;
        subplot(numel(b), 2, 2*i-1);
        plot(t, time_course(:,1), 'r', t, time_course(:,2)); % 横軸:時間、縦軸:濃度でプロット
        axis([350 time(end) 0 3]);
        
        %% figure 2 を描く
        hold on;
        subplot(numel(b),2,2:2:12);
        axis([0.0 2.5 0 3]);
        %plot_nullcline(param);
        ph(i)=plot(time_course(:,1), time_course(:,2), 'color', col(floor(56/numel(b)*i),:)); % 横軸:ADP、縦軸:F6Pでプロット
        leg1(i).name = ['b = ', num2str(b(i))];
%         waitforbuttonpress
        
    end
    
    %% figure 2 の凡例
    legend([ph],leg1.name);
    xlabel('ADP');
    ylabel('F6P');
end

function dsdt = ODE(t, s, param)
    x = s(1); % ADP
    y = s(2); % F6P
    a = param(1);
    b = param(2);
    
    v1 = 　　　　　　　　　　　　　; % 空欄を埋めよ
    v2 = 　　　　　　　　　　　　　; % 空欄を埋めよ
    v3 = 　　　　　　　　　　　　　; % 空欄を埋めよ
    
    dsdt(1, :) = 　　　　　　　　　　　　　; % 空欄を埋めよ % ADPの時間変化
    dsdt(2, :) = 　　　　　　　　　　　　　; % 空欄を埋めよ % F6Pの時間変化
end

function plot_nullcline(param)
    a = 　　　　　　　　　　　　　; % 空欄を埋めよ
    b = 　　　　　　　　　　　　　; % 空欄を埋めよ
    x = 0:0.05:2.5;
    y1 = 　　　　　　　　　　　　　; % 空欄を埋めよ
    y2 = 　　　　　　　　　　　　　; % 空欄を埋めよ
    plot(x, y1, 'k--', x, y2, 'k--');
end