%
% backward Euler on the fork
%
function befork(dt,Tfin,t1,t2,I0,pinc)

a = 1e-3*[1 1/2 2];
ell = [1 1 2];
dx = .025;
N = ell/dx;
A3 = 2*pi*a(3)*dx;
As = 4*pi*1e-6;
rho = A3/As;
R2 = 0.034;
gL = 0.3;
Cm = 1;
tau = Cm/gL;
lam = a/(2*R2*gL)/dx^2;   % lammda^2
r = a/a(3);
Hd = [2*lam(1)*ones(1,N(1)) 2*lam(2)*ones(1,N(2)) 2*lam(3)*ones(1,N(3)+1)];
Hd(1) = lam(1);
Hd(N(1)+1) = lam(2);
Hd(N(1)+N(2)+1) =  lam*r';
Hd(end) = rho*lam(3);
Hlen = length(Hd);

Hu = [-lam(1)*ones(1,N(1)-1) 0 -lam(2)*ones(1,N(2)) -lam(3)*ones(1,N(3))];
Hl = [-lam(1)*ones(1,N(1)-1) 0 -lam(2)*ones(1,N(2)-1) -r(2)*lam(2) -lam(3)*ones(1,N(3))];
Hl(end) = rho*Hl(end);

H = spdiags( [[Hl 0]' Hd' [0 Hu]'], -1:1, Hlen, Hlen);

H(N(1)+N(2)+1,N(1)) = -r(1)*lam(1);
H(N(1),N(1)+N(2)+1) = -lam(1);

I = speye(Hlen);

Bb = I+(I+H)*(dt/tau); 
[L,U] = lu(Bb); 

x3 = 0:dx:ell(3);
x1 = ell(3):dx:ell(3)+ell(1);
x2 = ell(3):dx:ell(3)+ell(2);

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

t = 0;
tcnt = 0;
x = linspace(0,ell(3)+max(ell(1:2)),Hlen);
plot3(x,t*ones(Hlen,1),v)
hold on

stimloc = N(1)+N(2)+5;

while (t <= Tfin)

      t = t + dt;

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

      rhs = v;
      rhs(Hlen) = rhs(Hlen) + dt*Istim/(Cm*As);
%      rhs(stimloc) = rhs(stimloc) + dt*Istim/(Cm*As);

      v = U \ (L \ rhs);

      tcnt = tcnt + 1;

      if mod(tcnt,pinc) == 0
         v1 = fliplr(v(1:N(1))');
         v2 = fliplr(v(N(1)+1:N(1)+N(2))');
         v3 = fliplr(v(N(1)+N(2)+1:end)');
         plot3(x3,t*ones(size(x3)),v3,'k')
         plot3(x2,t*ones(size(x2)),[v3(end) v2],'b')
         plot3(x1,t*ones(size(x1)),[v3(end) v1],'r')
      end

end

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