Dynamical system with bistability

deq = {u'[t] == u[t] - u[t]^3 - v[t],v'[t]==e(u[t] - a1 v[t] - a0),
u[0]==u0,v[0]==v0}

para = {a0->-0.1,a1->2,e->0.05}

deqN = deq/.para

tEnd = 100;
sol = NDSolve[deqN/.{u0->1,v0->-0.4},{u[t],v[t]},{t,0,tEnd}]

Plot[ Evaluate[{u[t],v[t]}/.sol],{t,0,tEnd},PlotRange->All,AxesLabel->{"t","u, v"}]

ginit1 = ParametricPlot[ Evaluate[{u[t],v[t]}/.sol],{t,0,tEnd},AxesLabel->{"u","v"},PlotRange->All,PlotStyle->{RGBColor[1,0,0]}]

tEnd = 100;
sol = NDSolve[deqN/.{u0->-1,v0->-0.4},{u[t],v[t]},{t,0,tEnd}]

Plot[ Evaluate[{u[t],v[t]}/.sol],{t,0,tEnd},AxesLabel->{"t","u, v"}]

ginit2 = ParametricPlot[ Evaluate[{u[t],v[t]}/.sol],{t,0,tEnd},AxesLabel->{"u","v"},PlotRange->All,PlotStyle->RGBColor[0,0,1]]

v1=  u - u^3/.para
v2 = (u-a0)/a1/.para

gphase = Plot[{v1,v2},{u,-2,2},PlotRange->{{-2,2},{-1,1}},PlotStyle->Thickness[0.005]]

Show[ginit1,ginit2,gphase,PlotRange->{{-2,2},{-1,1}}]

vec={u-u^3-v,e(u-a1 v-a0)}

vecN=vec/.para

gVecField=VectorPlot[vecN,{u,-1.5,1.5},{v,-1,1},VectorPoints->30]

Show[gphase,gVecField,PlotRange->{{-1.5,1.5},{-1,1}}]

Show[gphase,gVecField,ginit1,ginit2]

NSolve[Thread[vecN=={0,0}],{u,v}]

?VectorPlot