%
%  bepfib.m
%
%  Steve Cox, 1/20/03
%
% solve the passive fiber problem via backward euler
%
% usage:        bepfib(a,ell,N,dt,s,t1,t2,T,Iapp,pinc)
%
%               a = fiber radius (cm)
%               ell = fiber length (cm)
%               N = number of compartments
%               dt = timestep (ms)
%               s = current pulse compartment number
%               t1 = start of current pulse (ms)
%               t2 = end of current pulse (ms)
%               T = stopping time (ms)
%               I0 = size of current pulse (micro amps)
%               pinc = number of time steps between plots
%
% example:     bepfib(.001,1,100,.03,50,1,2,10,.01,10)
%

function bepfib(a,ell,N,dt,s,t1,t2,T,I0,pinc)

C_m = 1;		% micro F / cm^2
g_L = 0.3;     		% mS / cm^2
R_2 = 0.034;		% k Ohm cm
G_i = a/2/R_2;		% scaled axial conductance
dx = ell/N;		% patch length
A = 2*pi*a*dx;		% patch surface area
x = dx/2:dx:ell-dx/2;	% vector of patch midpoints

v = zeros(N,1);		% initial conditions
rhs = zeros(N,1);

S = (2*eye(N)-diag(ones(N-1,1),1)-diag(ones(N-1,1),-1))/dx/dx;
S(1,1) = 1/dx/dx;
S(N,N) = 1/dx/dx;

tau = C_m/g_L;
lambda = sqrt(a/(2*R_2*g_L));

B = (eye(N)+lambda^2*S)/tau;

Bb = eye(N) + dt*B;

t = 0;
tcnt = 0;
plot3(x,t*ones(N,1),v)
hold on

while (t <= T)

      t = t + dt;

      Istim = I0*(t>t1)*(t<t2);

      rhs = v;
      rhs(s) = rhs(s) + dt*Istim/A;
      v = Bb \ rhs;

      tcnt = tcnt + 1;

      if mod(tcnt,pinc) == 0
         plot3(x,t*ones(N,1),v)
      end

end

xlabel('x (cm)','fontsize',16)
ylabel('t (ms)','fontsize',16)
zlabel('v (mV)','fontsize',16)

hold off