%**************************************************************************
%*     Problem 5.2                                                        *
%*     filename: problem05_2.m                                            *
%*     Programed by Ryoichi Kawai for Computational Physics Course        *
%*     Revised on 012/07/2016.                                            *
%**************************************************************************
clear all;

% Get the value of I
% Somewhere between I=0.3241 and 0.3242, trhere is a bifurcation.
I=input('Enter the val;ue of I =');

% fixed parameter
a=0.08;
b=0.7;
c=0.8;

% right hand side of the odes
fv=@(v,w) v-v^3/3-w+I;
fw=@(v,w) a*(v+b-c*w);

%initial conditions
v(1)=1;
w(1)=1;
t(1)=0;

% duration 
tmax=300;

% number of integration steps
N=1000;

% step size
h=tmax/N;

for n=1:N-1
    % 4th-order Runge-Kutta
    kv1=fv(v(n),w(n));
    kw1=fw(v(n),w(n));
    
    v_mid = v(n)+kv1*h/2;
    w_mid = w(n)+kw1*h/2;
    kv2 = fv(v_mid,w_mid);
    kw2 = fw(v_mid,w_mid);
      
    v_mid = v(n)+kv2*h/2;
    w_mid = w(n)+kw2*h/2;
    kv3 = fv(v_mid,w_mid);
    kw3 = fw(v_mid,w_mid);

    v_end = v(n)+kv3*h;
    w_end = w(n)+kw3*h;
    kv4 = fv(v_end,w_end);
    kw4 = fw(v_end,w_end);
    
    v(n+1)=v(n)+(kv1+2*(kv2+kv3)+kv4)*h/6;
    w(n+1)=w(n)+(kw1+2*(kw2+kw3)+kw4)*h/6;
    t(n+1)=t(1)+n*h;
end

% plot individual trajectories
subplot(1,2,1);
p=plot(t,v,t,w);
xlabel('t');
ylabel('concentration');
set(p(1),'Color','blue','Linewidth',2);
set(p(2),'Color','red','Linewidth',2);
legend(p,{'v','w'});
legend(p,'Location','SouthEast');

% plot phase trajectory
subplot(1,2,2);
q=plot(v,w);
xlabel('v');
ylabel('w');
set(q(1),'color','blue');