Selkov-Oszillator

f[1,t_]=v0-k1*X[t]*Y[t]^2;
f[2,t_]=k1*X[t]*Y[t]^2-k2*Y[t];

(* stationärer Zustand *)
stat=Solve[{f[1,t]==0,f[2,t]==0},{X[t],Y[t]}]

k1=1;

(* Jacobi-Matrix *)
Jx=Table[D[f[i,t],X[t]]/.stat[[1]],{i,2}];
Jy=Table[D[f[i,t],Y[t]]/.stat[[1]],{i,2}];
J={Jx,Jy}

d=Det[J]
sp=Tr[J]
EW=Eigenvalues[J]

stab=Solve[sp==0,v0]

p1=Plot[Evaluate[v0/.stab[[2]]],{k2,0,10},PlotStyle->RGBColor[1,0,0]]

knfo=Solve[sp^2==4*d,v0]

p2=Plot[Evaluate[{v0/.knfo[[2]],v0/.knfo[[4]]}],{k2,0,10}]

Show[p1,p2,AspectRatio->1,Frame->True,FrameLabel->{k2,v0}]

k2=10;v0=80;
stat

int1=NDSolve[{X'[t]==f[1,t],Y'[t]==f[2,t],X[0]==1,Y[0]==8},{X[t],Y[t]},{t,0,10}];

Plot[Evaluate[X[t]/.int1],{t,0,10},PlotRange->All]
Plot[Evaluate[Y[t]/.int1],{t,0,10},PlotRange->All]

k2=10;v0=40;
stat

int2=NDSolve[{X'[t]==f[1,t],Y'[t]==f[2,t],X[0]==2,Y[0]==4},{X[t],Y[t]},{t,0,10}];

Plot[Evaluate[X[t]/.int2],{t,0,10},PlotRange->All]
Plot[Evaluate[Y[t]/.int2],{t,0,10},PlotRange->All]

ParametricPlot[Evaluate[{X[t],Y[t]}/.int2],{t,0,10}]

k2=10;v0=30;
stat

int3=NDSolve[{X'[t]==f[1,t],Y'[t]==f[2,t],X[0]==4,Y[0]==2},{X[t],Y[t]},{t,0,10}];

Plot[Evaluate[X[t]/.int3],{t,0,10},PlotRange->All]
Plot[Evaluate[Y[t]/.int3],{t,0,10},PlotRange->All]

ParametricPlot[Evaluate[{X[t],Y[t]}/.int3],{t,0,10}]

Nullcline1=Solve[{f[1,t]==0},{Y[t]}]
Nullcline2=Solve[{f[2,t]==0},{Y[t]}]

Plot[Evaluate[{((Sqrt[30])/(Sqrt[X])),(10/X),0}],{X,1,6}]

FinalPlot1 = ParametricPlot[Evaluate[{X[t],Y[t]}/.int3],{t,0,10}];

FinalPlot2 = Plot[Evaluate[{((Sqrt[30])/(Sqrt[X])),(10/X),0}],{X,1,6}];

Show[FinalPlot1,FinalPlot2]