Vt = [1:N]/N;

Vn= (-ga*(Vt-eps2).*(Vt-alf).*(Vt-1.)+Istim)./Vt ;
 Wn= - gs*(Vt-eps2).*(Vt-alf-1.) ;

wrange = [-.1 1.1 -.5 3];
