clear all
close all

global g v n a

g=2;
v=0.0285;
a=1;
n=0.1;


ci=[.3 .3]; % Condiciones iniciales

[t,y]=ode23(@selkov,[0 500],ci); % Integrador

su=y(:,1);
sd=y(:,2);

%%% Grafico de las nuclinas %%%
vv=0.08:.01:1;

su1=v./(1-v).*(1+vv.^g)./vv.^g;
su2=(1+vv.^g)./(vv.^(g-1).*(a/n-vv));

%%%%%%%%% calculo del pf %%%%%%%%%%%
sdf = a/n*v
suf = v*(1+sdf^g)/((1-v)*sdf^g)
%%%%%%%%%%%5

hold on
plot(suf,sdf,'*')
plot(su1,vv,'r',su2,vv,'b')




axis([min(su)*.9, max(su)*1.1, min(sd)*.9, max(sd)*1.1])
plot (su,sd,'g') % Grafico del diagrama de fases
xlabel('\sigma_1')
ylabel('\sigma_2')
hold off

figure
plot(t,su,'r',t,sd,'b') % Grafico de la evolucion temporal

xlabel('t')
ylabel('\sigma_1 y \sigma_2')