function F=driftvHR(q,v,stim)

global I C gbarNa ENa gbarK EK gbarL EL gA A

Iapp=I + stim;

F=(Iapp - gbarNa*minf(v)^3*(-3*(q-A*binf(v))+.85)*(v-ENa) - gbarK*q*(v-EK) -gbarL*(v-EL) )/C   ;