subroutine integrate(v,am,an,ah)
      implicit double precision (a-h,o-z)
      parameter(imax=500)
      dimension v(imax),am(imax),an(imax),ah(imax)
      dimension d(imax),b(imax),c(imax),r(imax)
      dimension chan(imax),pump(imax),pm(imax)
      common/initp/ dt,dx,v0,am0,an0,ah0
      common/params/ cm,gna,gk,gl,ena,ek,el,a,rho

      eki1 = a*dt/(dx*dx*2*rho*cm)
      eki2 = dt/cm
      do 10 i=2,imax-1
         alpham = (v(i)+45.0)/(10.0*(1.0-dexp((v(i)+45.0)/(-10.0)))) 
         betam = 4.0*dexp((v(i)+70.0)/(-18.0))
         alphan = 0.1*(v(i)+60.0)/(10.0*(1.0-dexp((v(i)+60.0)/(-10.0))))
         betan = 0.125*dexp((v(i)+70.0)/(-80.0))
         alphah = 0.07*dexp((v(i)+70.0)/(-20.0))
         betah = 1/(1+dexp((v(i)+40.0)/(-10.0)))
         taum = 1/(alpham+betam)
         aminf = alpham*taum
         taun = 1/(alphan+betan)
         aninf = alphan*taun
         tauh = 1/(alphah+betah)
         ahinf = alphah*tauh
         g = gna*(am(i)**3)*ah(i)+gk*(an(i)**4)+gl
         e = (gna*(am(i)**3)*ah(i)*ena+gk*(an(i)**4)*ek+gl*el)/g
         eki = g*eki2
         r(i) = v(i) + eki*e
         d(i) = -eki1 
         b(i) = 1. + eki + 2.*eki1
         c(i) = -eki1
         am(i) = (taum*am(i) + dt*aminf)/(taum + dt)
         an(i) = (taun*an(i) + dt*aninf)/(taun + dt)
         ah(i) = (tauh*ah(i) + dt*ahinf)/(tauh + dt)
   10 continue
      i=1
         alpham = (v(i)+45.0)/(10.0*(1.0-dexp((v(i)+45.0)/(-10.0))))
         betam = 4.0*dexp((v(i)+70.0)/(-18.0))
         alphan = 0.1*(v(i)+60.0)/(10.0*(1.0-dexp((v(i)+60.0)/(-10.0))))
         betan = 0.125*dexp((v(i)+70.0)/(-80.0))
         alphah = 0.07*dexp((v(i)+70.0)/(-20.0))
         betah = 1/(1+dexp((v(i)+40.0)/(-10.0)))
         taum = 1/(alpham+betam)
         aminf = alpham*taum
         taun = 1/(alphan+betan)
         aninf = alphan*taun
         tauh = 1/(alphah+betah)
         ahinf = alphah*tauh
         g = gna*(am(i)**3)*ah(i)+gk*(an(i)**4)+gl
         e = (gna*(am(i)**3)*ah(i)*ena+gk*(an(i)**4)*ek+gl*el)/g
         eki = g*eki2
         r(i) = v(i) + eki*e
         d(i) = 0.0
         b(i) = 1. + eki + eki1
         c(i) = -eki1
         am(i) = (taum*am(i) + dt*aminf)/(taum + dt)
         an(i) = (taun*an(i) + dt*aninf)/(taun + dt)
         ah(i) = (tauh*ah(i) + dt*ahinf)/(tauh + dt)
       i=imax
         alpham = (v(i)+45.0)/(10.0*(1.0-dexp((v(i)+45.0)/(-10.0))))
         betam = 4.0*dexp((v(i)+70.0)/(-18.0))
         alphan = 0.1*(v(i)+60.0)/(10.0*(1.0-dexp((v(i)+60.0)/(-10.0))))
         betan = 0.125*dexp((v(i)+70.0)/(-80.0))
         alphah = 0.07*dexp((v(i)+70.0)/(-20.0))
         betah = 1/(1+dexp((v(i)+40.0)/(-10.0)))
         taum = 1/(alpham+betam)
         aminf = alpham*taum
         taun = 1/(alphan+betan)
         aninf = alphan*taun
         tauh = 1/(alphah+betah)
         ahinf = alphah*tauh
         g = gna*(am(i)**3)*ah(i)+gk*(an(i)**4)+gl
         e = (gna*(am(i)**3)*ah(i)*ena+gk*(an(i)**4)*ek+gl*el)/g
         eki = g*eki2
         r(i) = v(i) + eki*e
         d(i) = -eki1
         b(i) = 1. + eki + eki1
         c(i) = 0.0
         am(i) = (taum*am(i) + dt*aminf)/(taum + dt)
         an(i) = (taun*an(i) + dt*aninf)/(taun + dt)
         ah(i) = (tauh*ah(i) + dt*ahinf)/(tauh + dt)
      call tridag(d,b,c,r,imax)
      do 100 i=1,imax
         v(i) = r(i)
  100 continue 
      return
      end