
%
%  cab1.m
%  Steve Cox, CAAM 335
%
%  Compute and plot the potentials in an N compartment
%  fiber resulting from a current stimulus at its left end
%

ell = 1;	% fiber length (cm)
a = 0.025;	% fiber radius (cm)
rho_i = 30;	% cytoplasmic resistivity (Ohm cm)
rho_m = 1000;	% membrane resistivity (Ohm cm^2)

N = 3;		% number of compartments

R_i = 10; %rho_i*(ell/N)/(pi*a*a);  	% axial resistance (Ohm)
R_m = 100; % rho_m/(2*pi*a*ell/N);	% lateral resistance (Ohm)

i_0 = .001;	% current stimulus (mA)

m = 2*N;	% number of edges
n = N+1;	% number of nodes

A = zeros(m,n);	% initialize adjacency matrix

for j=1:N,	% build A one compartment at a time

    A(2*j-1:2*j,j:j+1) = [-1 1;0 -1];

end

G = zeros(m,1);	% initialize and assemble G
G(1:2:m) = 1/R_i;
G(2:2:m) = 1/R_m;
G = diag(G);

S = A'*G*A;	% assemble main matrix

f = zeros(n,1);	% assemble the stimulus vector
f(1) = i_0;

x = S \ f;	% solve for the potentials

z = 0:ell/N:ell;	% nodal positions for plotting

plot(z,x,'x-')	% plot and label potential vs. distance
box off
xlabel('z  (cm)','fontsize',16)
ylabel('x  (mV)','fontsize',16)
