function [ output_args ] = selkov_phaseplane( input_args )
%
% Poblem
% Sel'kov model: 
% An example of biochemical oscillation and Hopf bifurcatino
%
    time = 0.001:1:50;
    s0 = [1.0 , 1.0];    % Initial values
    param = [0.06 , 0.6]; % Constants
    [t,time_course] = ode15s(@(t,s) ODE(t,s,param),time,s0); 
    
    plot_phase_plane(time_course,param(1),param(2));
end

function dsdt = ODE(t,s,param)
    ADP = s(1);
    F6P = s(2);
 
    a = param(1);
    b = param(2);
    
    v1 = ADP;
    v2 = a * F6P + ADP^2 * F6P;
    v3 = b
    
    dsdt(1,:) = - v1 + v2;
    dsdt(2,:) = v3 - v2;
end

function plot_phase_plane(time_course,a,b)
    figure;
    hold on;
    plot(        ,        );
    ADP = 0:0.1:3;
    F6P = ADP ./ ( a  + ADP.^2 ); % d ADP / dt == 0
    plot(ADP,F6P,'r');
    
    ADP = 0:0.1:3;
                     ;    % d F6P / dt == 0
    plot(      ,      ,      );
    hold off;
end
