function FHN
clear;
close;
time_end=50;
a=0.8;
b=0.9;
c=5;
z=input('z=');
%z=-0.5;
U=-2.5:0.01:2.5;
V=U.^3/3-U;
V1=(a-U)/b;
V2=U.^3/3-U-z;

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);
[T,Y] = ode45(@lveoms,[0 time_end],[1.3 -0.5],options);
plot(T,Y(:,1),'-b',T,Y(:,2),'-r','LineWidth',3);
axis([0 time_end -2.5 2.5]);
xlabel('t'); ylabel('u(blue), v(red)');
pause;
close;
[T,Y] = ode45(@lveoms,[0 time_end],[1.3 -0.5],options);
plot(U,V,'--g',U,V1,'--g',U,V2,'--r','LineWidth',3); 
hold on;
plot([-2.5 2.5],[0 0],'--k',[0 0],[-1.5 1.5],'--k','LineWidth',1);
axis([-2.5 2.5 -1.5 1.5]);
xlabel('u'); ylabel('v');
pause;
plot(Y(:,1),Y(:,2),'-b','LineWidth',3);
pause;
hold off;
close;

function dy=lveoms(t,y)
%  returns the coupled odes
    dy=zeros(2,1);
    dy(1)=c*(y(2)+y(1)-y(1)^3/3+z*0.5*(1+tanh((t-1.0)/0.1)));
    dy(2)=-(y(1)-a+b*y(2))/c;
end
end







 