minf = 0.5*(1.+tanh((V-v1)/v2));
 tmp = (V-v3)/v4; 
      winf = 0.5*(1.+tanh(tmp));
      r = cosh(tmp/2);
      current = gca*minf.*(V-vca)+gk*W.*(V-vk)+gl*(V-vl) -Iapp;
      Fv = -current/cm;
      Fw = phi*r.*(winf-W);
