% 2020-01-08

clearvars;
close all;

KSE = 3;

% Parameters

h = 2;
a = 3;
alpha = 4;
eps = 0.01;
lda = 0.008;
D = 0;

% Computation of the solution


T = 1000;
dt = 0.1;
t = 0:dt:T;

if KSE == 1

    v = zeros(1,length(t));
    w = zeros(1,length(t));

    wnse = randn(1,length(t));
    
    v(1) = 0;
    w(1) = 0.1;
    
    for j=1:length(t)-1
        k1v = -h*v(j)^3+a*v(j)^2-w(j);
        k1w = eps*(alpha*v(j)-lda-w(j));
        av = v(j) + k1v*dt;
        av = av + D*dt^0.5*wnse(j);
        aw = w(j) + k1w*dt;
        k2v = -h*av^3+a*av^2-aw;
        k2w = eps*(alpha*av-lda-aw);
        v(j+1) = v(j) + (k1v+k2v)*dt/2;
        v(j+1) = v(j+1) + D*dt^0.5*wnse(j);
        w(j+1) = w(j) + (k1w+k2w)*dt/2;
    end
    
    v(1) = v(end);
    w(1) = w(end);
    
    for j=1:length(t)-1
        k1v = -h*v(j)^3+a*v(j)^2-w(j);
        k1w = eps*(alpha*v(j)-lda-w(j));
        av = v(j) + k1v*dt;
        av = av + D*dt^0.5*wnse(j);
        aw = w(j) + k1w*dt;
        k2v = -h*av^3+a*av^2-aw;
        k2w = eps*(alpha*av-lda-aw);
        v(j+1) = v(j) + (k1v+k2v)*dt/2;
        v(j+1) = v(j+1) + D*dt^0.5*wnse(j);
        w(j+1) = w(j) + (k1w+k2w)*dt/2;
    end
    
    figure
    plot(t,v,'-b','linewidth',2);
    hold on
    axis([0 T -0.8 1.8])
    set(gca,'fontsize',20);
    xlabel('t  ');
    ylabel('v  ');
    
    vv = -2:0.001:2;
    vnlc = -h*vv.^3+a*vv.^2;
    wnlc = alpha*vv-lda;
    
    jin = floor(3*T/4/dt);
    jfin = floor(T/dt);
   
    
    figure
    vv=-2:0.001:2;
    plot(vv,vnlc,'-r','linewidth',2);
    hold on
    plot(vv,wnlc,'-g','linewidth',2);
    plot(v(jin:jfin),w(jin:jfin),'-b','linewidth',2);
    axis([-0.8 2 -0.1 1.4])
    set(gca,'fontsize',20);
    xlabel('v  ');
    ylabel('w');
    legend('v-nullcline','w-nullcline','trajectory');
    
    %axis([-0.2 0.2 -0.02 0.1])
elseif KSE == 2
    
    eta = 0.08;
    vth = 0.5;
    
    v = zeros(1,length(t));
    w = zeros(1,length(t));

    wnse = randn(1,length(t));
    
    v(1) = 0;
    w(1) = 0.1;
    z(1) = 0;
    
    lda = -0.1;
    
    for j=1:length(t)-1
        k1v = -h*v(j)^3+a*v(j)^2-w(j);
        k1w = eps*(alpha*v(j)-lda-z(j)-w(j));
        k1z = eps*eta;
        av = v(j) + k1v*dt;
        av = av + D*dt^0.5*wnse(j);
        aw = w(j) + k1w*dt;
        az = z(j) + k1z*dt;
        k2v = -h*av^3+a*av^2-aw;
        k2w = eps*(alpha*av-lda-z(j)-aw);
        k2z = eps*eta;
        v(j+1) = v(j) + (k1v+k2v)*dt/2;
        v(j+1) = v(j+1) + D*dt^0.5*wnse(j);
        w(j+1) = w(j) + (k1w+k2w)*dt/2;
        z(j+1) = z(j) + (k1z+k2z)*dt/2;
        if v(j+1) > vth 
            z(j+1) = 0;
        end
    end
    
    v(1) = min(v);
    w(1) = max(w);
    z(1) = 0;
   
    
    for j=1:length(t)-1
        k1v = -h*v(j)^3+a*v(j)^2-w(j);
        k1w = eps*(alpha*v(j)-lda-z(j)-w(j));
        k1z = eps*eta;
        av = v(j) + k1v*dt;
        av = av + D*dt^0.5*wnse(j);
        aw = w(j) + k1w*dt;
        az = z(j) + k1z*dt;
        k2v = -h*av^3+a*av^2-aw;
        k2w = eps*(alpha*av-lda-z(j)-aw);
        k2z = eps*eta;
        v(j+1) = v(j) + (k1v+k2v)*dt/2;
        v(j+1) = v(j+1) + D*dt^0.5*wnse(j);
        w(j+1) = w(j) + (k1w+k2w)*dt/2;
        z(j+1) = z(j) + (k1z+k2z)*dt/2;
        if v(j+1) > vth 
            z(j+1) = 0;
        end
    end
    
    
    
    
    
    figure
    plot(t,v,'-b','linewidth',2);
    hold on
    axis([0 T -0.8 1.8])
    set(gca,'fontsize',20);
    xlabel('t  ');
    ylabel('v  ');
    
    vv = -2:0.001:2;
    vnlc = -h*vv.^3+a*vv.^2;
    wnlc = alpha*vv-lda;
    
    jin = floor(3*T/4/dt);
    jfin = floor(T/dt);
   
    
    figure
    vv=-2:0.001:2;
    plot(vv,vnlc,'r','linewidth',2);
    hold on
    plot(vv,wnlc,'g','linewidth',2);
    plot(v(jin:jfin),w(jin:jfin),'b','linewidth',2);
    axis([-0.8 2.2 -0.1 1.5])
    set(gca,'fontsize',20);
    xlabel('v  ');
    ylabel('w');
    legend('V-nullcline','w-nullcline','trajectory');
elseif KSE == 3
    
    eta = 0.05;
    beta = alpha;
    sgma = -0.9;
    vth = 0.5;
    
    v = zeros(1,length(t));
    w = zeros(1,length(t));

    wnse = randn(1,length(t));
    
    v(1) = 0;
    w(1) = 0.1;
    z(1) = 0;
    
    lda = -0.02;
    
    for j=1:length(t)-1
        k1v = -h*v(j)^3+a*v(j)^2-w(j);
        k1w = eps*(alpha*v(j)-lda-z(j)-w(j));
        k1z = eps*eta*(beta*v(j)-sgma-z(j));
        av = v(j) + k1v*dt;
        av = av + D*dt^0.5*wnse(j);
        aw = w(j) + k1w*dt;
        az = z(j) + k1z*dt;
        k2v = -h*av^3+a*av^2-aw;
        k2w = eps*(alpha*av-lda-z(j)-aw);
        k2z = eps*eta*(beta*av-az);
        v(j+1) = v(j) + (k1v+k2v)*dt/2;
        v(j+1) = v(j+1) + D*dt^0.5*wnse(j);
        w(j+1) = w(j) + (k1w+k2w)*dt/2;
        z(j+1) = z(j) + (k1z+k2z)*dt/2;
        if v(j+1) > vth 
            z(j+1) = 0;
        end
    end
    
    v(1) = min(v);
    w(1) = max(w);
    z(1) = 0;
   
    
    for j=1:length(t)-1
        k1v = -h*v(j)^3+a*v(j)^2-w(j);
        k1w = eps*(alpha*v(j)-lda-z(j)-w(j));
        k1z = eps*eta*(beta*v(j)-sgma-z(j));
        av = v(j) + k1v*dt;
        av = av + D*dt^0.5*wnse(j);
        aw = w(j) + k1w*dt;
        az = z(j) + k1z*dt;
        k2v = -h*av^3+a*av^2-aw;
        k2w = eps*(alpha*av-lda-z(j)-aw);
        k2z = eps*eta*(beta*av-sgma-az);
        v(j+1) = v(j) + (k1v+k2v)*dt/2;
        v(j+1) = v(j+1) + D*dt^0.5*wnse(j);
        w(j+1) = w(j) + (k1w+k2w)*dt/2;
        z(j+1) = z(j) + (k1z+k2z)*dt/2;
        if v(j+1) > vth 
            z(j+1) = 0;
        end
    end
    
  
    
    
    figure
    plot(t,v,'linewidth',2);
    hold on
    axis([0 T -0.8 1.8]);
    set(gca,'fontsize',20);
    xlabel('t  ');
    ylabel('v  ');
    
    vv = -2:0.001:2;
    vnlc = -h*vv.^3+a*vv.^2;
    wnlc = alpha*vv-lda;
    
    jin = floor(3*T/4/dt);
    jfin = floor(T/dt);
   
    
    figure
    vv=-2:0.001:2;
    plot(vv,vnlc,'r','linewidth',2);
    hold on
    plot(vv,wnlc,'g','linewidth',2);
    plot(v(jin:jfin),w(jin:jfin),'b','linewidth',2);
    axis([-0.8 2.2 -0.1 1.5])
    set(gca,'fontsize',20);
    xlabel('v  ');
    ylabel('w');
    legend('V-nullcline','w-nullcline','trajectory');

end
        
        
    
    