An extended Hodgkin-Huxley model that doesn't want to work.

Hi Berend

This is the code. Pretty much just changed to what you suggested which is
CaI=1, removing the unnecessary variables and using deSolve:

rm(list = ls())
library(deSolve)
## Hodkin-Huxley model
HH_soma <-  function(times, init, parms) {
 with(as.list(c(init, parms)),{

   # Na only used in Axon
   #Naminf <-1/(1+exp(-(v+24.7)/5.29));
   #Nataum <- function(v) 1.32 - (1.26/(1+exp(-(v+120)/25)));
   #Nahinf <-1/(1+exp((v+489)/5.18));
   #Natauh <-(0.67/(1+exp(-(v+62.9)/10))) * (1.5+(1/(1+exp((v+34.9)/36))));

   #PD
   # mca10
   CaTminf <- function(v) 1/(1+exp(-(v+25)/7.2));
   # hca10
   CaThinf <- function(v) 1/(1+exp(v+36)/7);
   # taumca1
   CaTtaum <- function(v) 55- (49.5/(1+exp(-v+58)/17));
   # tauhca1
   CaTtauh <- function(v) 350 - (300/(1+exp(-v+50)/16.9));

   #mca20
   CaSminf <- function(v) 1/(1+exp(-(v+22)/8.5));
   #taumca2
   CaStaum <- function(v) 16-(13.1/(1+exp(-(v+25.1)/26.4)));


   # mna0
   Napminf <- function(v) 1/(1+exp(-(v+26.8)/8.2));
   # hna0
   Naphinf <- function(v) 1/(1+exp(-(v+48.5)/5.18));

   taumna <- function(v) 19.8-(10.7/(1+exp(-(v+26.5)/8.6)));
   tauhna <- function(v) 666-(379/(1+exp(-(v+33.6)/11.7)));

   # mh0
   hminf <- function(v)  1/(1+exp(v+70)/6);
   # taumh
   htaum <- function(v)  272+(1499/(1+exp(-(v+42.2)/8.73)));

   Kminf <- function(v)  1/(1+exp(-(v+14.2)/11.8));
   Ktaum <- function(v)  7.2-(6.4/(1+exp(-(v+28.3)/19.2)));

   # Reversal potential of intracellular calcium concentration
   # Nernst Equation using  extracellular concentration of Ca = 13mM
   # eca
   ECa <- function(Ca2) 12.2396*log(13000/(Ca2));
   #ECa <- function(CaI) 12.2396*log(13000/(CaI));


   #Sum of all the Ca
   # function(v) CaTminf(v) + CaSminf(v);
   CaI <- gCaT*CaTm^3*CaTh*(v-ECa(CaI)) + gCaS*CaSm^3*(v-ECa(CaI))

   #AB
   #dCa2 <- (((-F*Caminf(v))-Caminf(v) + C0)/TauCa)
   dCa2 <- (((-F*CaI) - Ca2 + C0)/TauCa)

   # mk20
   KCaminf <- function(v, Ca2) (Ca2/(Ca2+30))*(1/(1+exp(-(v+51)/8)));
   # taumk
   KCataum <- function(v) 90.3 - ((75.09/(1+exp(-(v+46)/22.7))));

   #AB
   Aminf <- function(v) 1/(1+exp(-(v+27)/8.7));
   Ahinf <- function(v) 1/(1+exp((v+56.9)/4.9));
   Ataum <- function(v) 11.6-(10.4/(1+exp(-(v+32.9)/15.2)));
   Atauh <- function(v) 38.6-(29.2*(1+exp(-(v+38.9)/26.5)));

   #proc
   #mp0
   procminf <- function(v) 1/(1+exp((v+56.9)/4));
   #taump
   proctaum <- function(v) 0.5;


    dv <- (-1*(I
              + CaI
              + gNap*Napm^3*Naph*(v-ENap)
              + gh*hm*(v-Eh)
              + gK*Km^4*(v-EK)
              + gKCa * KCam^4*(v-EKCa)
              + gA*Am^4*Ah*(v-EA)
              + gL*(v-EL))
          / C);


   dCaTm <- (CaTminf(v) - CaTm)/CaTtaum(v);
   dCaTh <- (CaThinf(v) - CaTh)/CaTtauh(v);

   dCaSm <- (CaSminf(v) - CaSm)/CaStaum(v);

   dNapm <- (Napminf(v) - Napm)/taumna(v);
   dNaph <- (Napminf(v) - Naph)/tauhna(v);

   dhm <- (hminf(v) - hm)/htaum(v);

   dKm <- (Kminf(v) - Km)/Ktaum(v);

   dKCam <- (KCaminf(v, Ca2) - KCam)/KCataum(v);

   dAm <- (Aminf(v) - Am)/Ataum(v);
   dAh <- (Ahinf(v) - Ah)/Atauh(v);


   list(c(dv,
          dCaTm, dCaTh,
          dCaSm,
          dNapm, dNaph,
          dhm,
          dKm,
          dKCam,
          dAm, dAh))
 })
}
parms = c(gCaT=22.5, gCaS=60, gNap=4.38, gh=0.219,
         gK=1576.8, gKCa=251.85, gA=39.42, gL=0.105,
         ENap=50, Ca2=0.52,
         Eh=-20, EK=-80, EL=-55, EKCa=-80, EA=-80,
         C=1/12, I=6.5, F=0.418, TauCa=303, C0=0.5, CaI=1);
time = seq(from=0, to=1000, by=0.1);
init = c(v=-65, CaTm=0.52 , CaTh=0.52, CaSm=0.52,
        Napm=0.52, Naph=0.52, hm=0.52, Km=0.52,
        KCam=0.52, Am=0.52, Ah=0.52);


out<-ode(y=init, times=time, func=HH_soma, parms=parms);
plot(ode)
o<-data.frame(out);
plot(o$time, o$v, type='l');