%%
%% Bifurcation in parameter eps = C/L --- measuring capacitor and inductor's energy generation capability. Small eps, large potiential energy.
%% eps = C/L --- measuring capacitor and inductor's 
%%                energy generation capability. 
%%                Small eps, large potiential energy.
%%
%%

clear all
global a b c  eps R E I_in

a=2; b=-2; c=1.1; 
%
% for the nonlinear resistor's IV-characteristic: d*V.*(a*V.^2+b*V+c);
%

R=0.5; E=.4; I_in=-.6; 
%%
%% R--resistance, E--battery 
%% eps = C/L with L--inductance, C--capacitance, 
%% 

% I_in=0; % stop the oscillation. 
% E = 0; % stop the oscillation.
eps = .001; 
% eps > 1 % stop the oscillation.

%%
%% Bifurcation diagram 
%%
xinit=[0.8243; 0.2325];
epsa=.001;
epsb=1.2;
K=40;
epsinc=(epsb-epsa)/K;
Vextrema=[];
Iextrema=[];
for k=0:K
   eps = epsa+k*epsinc;
   [t,x]=ode15s('FHN',[0,20],xinit);
   xinit=x(end,:); % Adaptive initials for different eps: current eps' orbital end is the next eps' orbital initial.
   N=find(t>2); % Throw away transient state and only keep the last t=2 asymptotic state.
   y=x(N(1):end,:);
   Vextrema=[Vextrema;[min(y(:,1)),max(y(:,1))]]; % max and min of V_C
   Iextrema=[Iextrema;[min(y(:,2)),max(y(:,2))]]; % max and min of I_L
end
epsbif=epsa:epsinc:epsb;
hold off
plot(epsbif',Vextrema(:,1),'r-',epsbif',Vextrema(:,2),'r-','linewidth',3)
hold on; plot([epsb epsb+.3], Vextrema(end,:),'r-','linewidth',3) % Extend the equilibrium branch
xlabel('\it\epsilon','fontsize', 24)
ylabel('Range of \it{V_C}','fontsize', 18)




