NumbEqus:=9:
NumbBCs:=12:  # + 3 unknown params

#...The Equations that we need to solve, formulated as a first order system
#-----------------------------3---------------------------
 F[1] := diff(E3(z),z)     = l*(Eprim3(z)):
 F[2] := diff(Eprim3(z),z) = l*(-c*Eprim3(z)):
 F[3] := diff(h3(z),z)     = l*(-1/c/th(E3(z))*(1-h3(z))):
#-----------------------------2---------------------------
 F[4] := diff(E2(z),z)     = (-l/p)*Eprim2(z):
 F[5] := diff(Eprim2(z),z) = (-l/p)*(-c*Eprim2(z)):
 F[6] := diff(h2(z),z)     = (-l/p)*(-1/c/th(E2(z))*(-h2(z))):
#-----------------------------1---------------------------
 F[7] := diff(E1(z),z)     = (-l)*Eprim1(z):
 F[8] := diff(Eprim1(z),z) = (-l)*(-c*Eprim1(z) -gNa*jalpha*(eNa-E1(z))*h1(z)):
 F[9] := diff(h1(z),z)     = (-l)*(-1/c/th(E1(z))*(-h1(z))):
#---------------------------------------------------------

#The boundary conditions in the format of the D02RAF are
 AL[1] := E3(0)      = Vh:
 AL[2] := E2(0)      = Vh:
 AL[3] := E1(0)      = Vm:
 AL[4] := Eprim3(0)  = Eprim2(0):
 AL[5] := h3(0)      = h2(0):
 AL[6] := Eprim1(0)  = Eprim2(1):
 AL[7] := h1(0)      = h2(1):
 AL[8] := Eprim3(1)  = -c*(E3(1)-Ealpha):
 AL[9] := E2(1)      = Vm:
 AL[10] := Eprim1(1) = (E1(1)-Eomega)/c/th(E1(1)):
 AL[11] := h3(1)     = 1:
 AL[12] := h1(1)     = -Eprim1(1)/gNa/jalpha/(eNa-Eomega)*(1/c/th(Eomega)+c):

#--- constants
eNa := 50.:
eK  := -85:
gNa := 4.:
gK  := 0.35:
Vm:=-33.75:
Vh:=-71.33:

l:=10:

BVPROBLEM:=[seq(F[i],i=1..NumbEqus), seq(AL[i],i=1..NumbBCs)]:

th:=proc(V)
 local ah, bh:
 ah := (0.1260*exp(-0.250*(V+77.))):
 bh := 1.7/(exp(-0.082*(V+22.5))+1.):
 1./(ah+bh);
end proc:
