% This subroutine does the LU decomposition for the
% orthogonal collocation method on finite elements.
% This version is for a single unknown, giving a block
% diagonal matrix with a single entry of overlap.
% This subroutine is in the Appendix
% of "Nonlinear Analysis in Chemical Engineering".
% It assumes diagonal pivoting is OK.
% The matrix Aocfe is stored as Aocfe(np,np,ne), where
% np is the number of collocation points in one element
% and ne is the number of elements.

function OCFElud(np,ne)
global Aocfe

n1 = np - 1;
for l=1:ne
	for k=1:n1
		k1=k+1;
		for i=k1:np
			s = Aocfe(i,k,l)/Aocfe(k,k,l);
			Aocfe(i,k,l) = s;
			for j=k1:np
				Aocfe(i,j,l) = Aocfe(i,j,l) - s*Aocfe(k,j,l);
			end
		end
	end
	if l<ne
		Aocfe(1,1,l+1) = Aocfe(np,np,l);
	end
end


% This subroutine does the fore and aft sweep for the
% orthogonal collocation method on finite elements.
% This version is for a single unknown, giving a block
% diagonal matrix with a single entry of overlap.
% This subroutine is in the Appendix
% of "Nonlinear Analysis in Chemical Engineering".
% It assumes diagonal pivoting is OK.
% The matrix Aocfe is stored as Aocfe(np,np,ne), where
% np is the number of collocation points in one element
% ne is the number of elements.
% nt is the total number of unknowns, nt = (np-1)*ne + 1
% To use OCFEfas, you must first run OCFElud.
% The right-hand side is in Bocfe(nt)

function OCFEfas(np,ne,nt)
global Aocfe Bocfe

np1=np;
% forward sweep
for l=1:ne
	for i=2:np1
		i2=(l-1)*(np-1)+i;
		s = 0.;
		i1=i-1;
		for j=1:i1
			j2=i2-i+j;
			s=s+Aocfe(i,j,l)*Bocfe(j2);
		end
		Bocfe(i2)=Bocfe(i2)-s;
	end
end
% back substitution
for l1=1:ne
	l=ne-l1+1;
	if (l==ne)
		Bocfe(nt)=Bocfe(nt)/Aocfe(np,np,ne);
	end
	n1 = np-1;
	for k=1:n1
		i=n1+1-k;
		i2=(l-1)*(np-1)+i;
		m=i+1;
		n2=n1+1;
		s=0.;
		for j=m:n2
			j2=(l-1)*(np-1)+j;
			s=s+Aocfe(i,j,l)*Bocfe(j2);
		end
		Bocfe(i2)=(Bocfe(i2) - s)/Aocfe(i,i,l);
	end
end




function lu_tri(n)
global afd bfd cfd

%     this subroutine does an lu decomposition of
%	  a tridiagonal system of 
%     equations of the type that often occur with the
%     finite difference method.
%	  After calling lu_tri one must call fas_tri to
%	  process the right-hand side.
%     input
%     a(n), b(n), c(n) in the following equation:
%     a(i)*y(i-1) + b(i)*y(i) + c(i)*y(i+1) = d(i)
%     n - the size of the system being solved.

%   output
%     afd(n), bfd(n), cfd(n) this is the lu decomposition.
%     the original afd, bfd, cfd are destroyed.

%  lower decomposition  tridiagonal matrix

for l=2:n
  s = afd(l)/bfd(l-1);
  bfd(l) = bfd(l)-s*cfd(l-1);
  afd(l) = s;
end


function fas_tri(n)
global afd bfd cfd dfd

%     this subroutine does the fore and aft
%	  sweep to solve a tridiagonal system of	   
%     equations of the type that often occur with the
%     finite difference method.
%	  One must call lu_tri one must call first.
%     input
%     afd(n), bfd(n), cfd(n) from the lu decomposition
% 	  solving the problem:
%     a(i)*y(i-1) + b(i)*y(i) + c(i)*y(i+1) = d(i)
%     n - the size of the system being solved.
%   output
%     dfd(n) is the solution

%  forward sweep  tridiagonal matrix
for l=2:n
   dfd(l) = dfd(l)-afd(l)*dfd(l-1);
end
%  back substitution
dfd(n) = dfd(n)/bfd(n);
for l=2:n
   k = n-l+1;
   dfd(k) = (dfd(k)-cfd(k)*dfd(k+1))/bfd(k);
end