with(DEtools): with(plots):

#Flow rate model
V0:=0;Qin:=1; Vb:=V0+Qin*t;
gstretch:=piecewise(t<1100,gsyn*sig(Vb),0);
sig:=t->1/(1+exp(0.01*(-t+400)));
#plot(gsyn*sig(t),t=0..1100);
#quit;

DE := 
diff(v(t), t)=-v(t)+Heaviside(v(t)-a)   -w(t) -gstretch*(v(t)-vsyn), 
diff(w(t), t) =eps*(v(t)+d-c*w(t));

a := .3; eps := .08; c := -0.3; d:=-1/5;
gsyn := -0.6; vsyn := 1/2;

solve({-ve+Heaviside(ve-a)-we, eps*(ve+d-c*we)}): assign(%);
IC1 := v(0) =ve, w(0) =we;


A1 :=plot(1/c*(v+d), v = -2 .. 2.1, color = blue):
#BB1 :=plot(-v+Heaviside(v-a), v = -2 .. 2.1, color = red):
BB1  :=plot(-v+Heaviside(v-a)-gsyn*(v-vsyn), v = -2 .. 2.1, color = red):

p1 := dsolve({DE, IC1}, {v(t),w(t)},numeric): 
#odeplot(p1, [t, w(t)], 0 .. 1000);
  odeplot(p1, [t, v(t)], 0 .. 1600,numpoints=4000,thickness=2);

#quit;

C1:=phaseportrait({DE},
               [v(t), w(t)],
               t=0..2000,
		 {[0,ve,we]},
                v =-1 .. 2.1,w = -0.5 ..0.9,
                arrows = SLIM, color=black,
                stepsize =.1,scene = [v,w]):
display(A1,BB1,C1);

