I = 0;
a = 1/4;
b = 1/3;
eps = 0.001;

function [f] = FHN(t,X)
    u = X(1)
    v = X(2)
    f(1) = u*(1-u)*(u-a)-v+I
    f(2) = eps*(u-b*v)
    f = [f(1);f(2)]
endfunction

dt = 0.01;
dx = 0.05;
Tf = 2000;
t = 0:dt:Tf;
x = -0.25:dx:1.2;

subplot(121)
for k=0.05:0.05:1.2
    u0 = [k;0];
    [u] = ode(u0,0,t,FHN)
    plot2d(u0(1),u0(2),-9,"000")
    plot2d(u(1,:),u(2,:),2)
    p = get("hdl");
    p.children.thickness = 2;
end

plot2d(x,x.*(1-x).*(x-a)+I,3)
p = get("hdl");
p.children.thickness = 3;
a=get("current_axes")
a.font_size=4; 
a.data_bounds=[-0.25,-0.02;1.2,0.12];
xtitle("Portrait de phase")

subplot(122)
plot2d(t,u(1,:),1)
plot2d(t,u(2,:),2)
p = get("hdl");
p.children.thickness = 3;
xtitle("Une trajectoire")
