function fitzhugh_nagumo_parameter_dependency()
    % selcovモデル
    % ODEを数値的に解く
    %% パラメータの指定
    s0 = [1, 0]; % 膜電位(v)の初期値と不活性化変数(w)の初期値
    param = [0.0 , 1, 0, 10]; % パラメータ a, b, Input, tau
    time = 0.01:0.05:300; %シミュレーションを行う時間
    Input = -0.4:0.1:0.4; % パラメータを振る
    col = jet; % 色指定に使用
    figure(1);
    clf;
    figure(2);
    clf;
    
    %% パラメータを振り分けて計算
    for i = 1:numel(Input);
        param(3) = Input(i); %% 振ったパラメータを代入
        
        %% 数値計算
        [t, time_course] = ode45(@(t, s) ODE(t, s, param), time, s0); % ODEを解く
        
        %% figure 1 を描く
        figure(1);
        hold on;
        subplot(numel(Input), 1, i);
        plot(t, time_course(:,1), 'r', t, time_course(:,2)); % 横軸:時間、縦軸:濃度でプロット
        axis([0 time(end) -2 2]);
        
        %% figure 2 を描く
        figure(2);
        hold on;
        axis([-2 2 -1.2 1.2]);
        plot_nullcline(param);
        ph(i)=plot(time_course(:,1), time_course(:,2), 'color', col(floor(53/numel(Input)*i),:)); % 横軸:時間、縦軸:濃度でプロット
        leg1(i).name = ['Input = ', num2str(Input(i), '%+6.2f')];
        
    end
    
    %% figure 2 の凡例
    figure(2);
    legend([ph],leg1.name,'Location','NorthEastOutside');
    xlabel('v');
    ylabel('w');
end

function dsdt = ODE(t, s, param)
    v = s(1); % 膜電位
    w = s(2); % 不活性化変数
    a = param(1);
    b = param(2);
    Input = param(3);
    tau = param(4);
    
    dsdt(1, :) = v - v.^3./3 - w + Input; % 膜電位の時間変化
    dsdt(2, :) = 1 ./ tau .* (v - a - b .* w); % 不活性化変数の時間変化
end


function plot_nullcline(param)
    a = param(1);
    b = param(2);
    Input = param(3);
    tau = param(4);
    v = -2.5:0.05:2.5;
    w1 = v - v.^3./3 + Input;
    w2 = 1./b .* (v - a);
    plot(v, w1, 'k', v, w2, 'k');
end