%
% fib3.m, Steve Cox, August 1, 2000
% code associated with chapter 6
%

A = [-1 0 0
     -1 0 0
     -1 1 0
     0 -1 0
     0 -1 0
     0 -1 1
     0 0 -1
     0 0 -1];

ell = .1;        % fiber length (cm)
a = 0.025;      % fiber radius (cm)
rho_i = 100;     % cytoplasmic resistivity (Ohm cm)
rho_m = 1000;   % membrane resistivity (Ohm cm^2)
Em = 0;      % resting potential mV
N = 2;         % number of compartments
c = .1;		% micro farad
R_i = rho_i*(ell/N)/(pi*a*a);   % axial resistance (Ohm)
R_m = rho_m/(2*pi*a*ell/N);     % lateral resistance (Ohm)

A_cb = 4*pi*(.05)^2;
R_cb = rho_m/A_cb;
Ccb = c*A_cb;
Cm = c*(2*pi*a*ell/N);

G = diag([0 1/R_cb 1/R_i 0 1/R_m 1/R_i 0 1/R_m]);

C = diag([Ccb 0 0 Cm 0 0 Cm 0]);

b = -Em*[0 1 0 0 1 0 0 1].';

B = -inv(A'*C*A)*(A'*G*A);

dt = 1;

S = (eye(3)/dt-B);

x = zeros(3,400);

t(1) = 0;

for j=2:400/dt,

    t(j) = (j-1)*dt;
    g = [t(j)^3*exp(-t(j)/6)/10000/Ccb 0 0]';
    x(:,j) = S \ ( x(:,j-1)/dt + g);

end

plot(t,x')
legend('x1','x2','x3')
xlabel('t (ms)','fontsize',16)
ylabel('x (mV)','fontsize',16)