% solve_CN_tridiagonal_system.m
% Procedure called by CNfd_Amer_put.m

% substitute boundary condition at j=-Nj into j=-Nj+1:

pmp(-Nj+1+joffset) = pm + pd;
pp(-Nj+1+joffset) = -pu*C(0+ioffset,-Nj+2+joffset)-(pm-2)*C(0+ioffset,-Nj+1+joffset) - pd*C(0+ioffset,-Nj+joffset) + pd*lambda_L;

% eliminate upper diagonal:

for j=-Nj+2:Nj-1
    pmp(j+joffset) = pm -pu*pd/pmp(j-1+joffset);
    pp(j+joffset) = -pu*C(0+ioffset,j+1+joffset) -(pm-2)*C(0+ioffset,j+joffset)-pd*C(0+ioffset,j-1+joffset) - pp(j-1+joffset)*pd/pmp(j-1+joffset);
end % for j

% use boundary condition at j=Nj and equation at j=Nj-1:

C(1+ioffset,Nj+joffset) = (pp(Nj-1+joffset) + pmp(Nj-1+joffset)*lambda_U)/(pu + pmp(Nj-1+joffset));
C(1+ioffset,Nj-1+joffset) = C(1+ioffset,Nj+joffset) - lambda_U;

% back-substitution:

for j=Nj-2:-1:-Nj+1
    C(1+ioffset,j+joffset) = ( pp(j+joffset) - pu*C(1+ioffset,j+1+joffset) )/pmp(j+joffset);
end % for j
% for j=-Nj, use equation (3.41) of text:
C(1+ioffset,-Nj+joffset) = C(1+ioffset,-Nj+1+joffset)-lambda_L;

% end of procedure solve_implicit_tridiagonal_system.m