function coll(nx,aaa,kon)
% coll - for orthogonal collocation with symmetric polynomials
% nx+1 = total number of collocation points
% aaa = geometry factor (1 = planar, 2 = cylindrical, 3 = spherical)
% kon = 0, will print matrices
% xxcoll(I) = location of collocation points, X(1) >0, X(N+1) = 1.0
% The collocation points for nx = 1 are for W = 1-x^2
% The collocation points for nx > 1 are for W = 1
% returned output is:
%	xxcoll(i), the collocation points
%	aacoll(i), the matrix for first derivative
%	bbcoll(i), the matrix for the Laplacian
%	qinvcoll(i), the matrix for the inverse of q

global xxcoll aacoll bbcoll qinvcoll wwcoll

if(nx<1) 
  nx=1;
  fprintf('\n ***************************************** \n')
  fprintf('The number of interior collocation points must be between \n')
  fprintf('1 and 6 (inclusive).  nx has been changed to %g\n\n',nx)
end
if(nx>6)
  nx=6;
  fprintf('\n ***************************************** \n')
  fprintf('The number of interior collocation points must be between \n')
  fprintf('1 and 6 (inclusive).  nx has been changed to %g\n\n',nx)
end

% now the code
n1 = nx+1;
%fprintf('The total number of collocation points (symmetric) is %g\n\n',n1)
%fprintf('The geometry factor is (1=planar, 2=cyl., 3=sph.) %g\n\n',aaa)
xxcoll(n1) = 1.0;

% find the collocation points
if (aaa==3)    %spherical
	if nx == 1
		xxcoll(1) = .6546536707;
	elseif nx==2
		xxcoll(1) = .5384693101;
		xxcoll(2) = .9061798459;
	elseif nx==3
		xxcoll(1) = .4058451514;
		xxcoll(2) = .7415311856;
		xxcoll(3) = .9491079123;
	elseif nx==4
		xxcoll(1) = .3242534234;
		xxcoll(2) = .6133714327;
		xxcoll(3) = .8360311073;
		xxcoll(4) = .9681602395;
	elseif nx==5
		xxcoll(1) = .2695431560;
		xxcoll(2) = .5190961292;
		xxcoll(3) = .7301520056;
		xxcoll(4) = .8870625998;
		xxcoll(5) = .9782286581;
	elseif nx==6
		xxcoll(1) = .2304583160;
		xxcoll(2) = .4484927510;
		xxcoll(3) = .6423493394;
		xxcoll(4) = .8015780907;
		xxcoll(5) = .9175983992;
		xxcoll(6) = .9841830547;
	end	
elseif (aaa==2)   %cylindrical
	if nx == 1
		xxcoll(1) = .5773502692;
	elseif nx==2
		xxcoll(1) = .4597008434;
		xxcoll(2) = .8880738340;
	elseif nx==3
		xxcoll(1) = .3357106870;
		xxcoll(2) = .7071067812;
		xxcoll(3) = .9419651451;
	elseif nx==4
		xxcoll(1) = .2634992300;
		xxcoll(2) = .5744645143;
		xxcoll(3) = .8185294874;
		xxcoll(4) = .9646596062;
	elseif nx==5
		xxcoll(1) = .2165873427;
		xxcoll(2) = .4803804169;
		xxcoll(3) = .7071067812;
		xxcoll(4) = .8770602346;
		xxcoll(5) = .9762632447;
	elseif nx==6
		xxcoll(1) = .1837532119;
		xxcoll(2) = .4115766111;
		xxcoll(3) = .6170011402;
		xxcoll(4) = .7869622564;
		xxcoll(5) = .9113751660;
		xxcoll(6) = .9829724091;
	end	
else   %planar
	if nx == 1
		xxcoll(1) = .4472135955;
	elseif nx==2
		xxcoll(1) = .3399810436;
		xxcoll(2) = .8611363116;
	elseif nx==3
		xxcoll(1) = .2386191861;
		xxcoll(2) = .6612093865;
		xxcoll(3) = .9324695142;
	elseif nx==4
		xxcoll(1) =.1834346425;
		xxcoll(2) =.5255324099; 
		xxcoll(3) =.7966664774;
		xxcoll(4) =.9602898565;
	elseif nx==5
		xxcoll(1) = .1488743390;
		xxcoll(2) = .4333953941;
		xxcoll(3) = .6794095683;
		xxcoll(4) = .8650633667;
		xxcoll(5) = .9739065285;
	elseif nx==6
		xxcoll(1) = .1252334085;
		xxcoll(2) = .3678314990;
		xxcoll(3) = .5873179543;
		xxcoll(4) = .7699026742;
		xxcoll(5) = .9041172564;
		xxcoll(6) = .9815606342;
	end
end

% form Q matrix
for i=1:n1
    z(i) = 1./(2.*i+aaa-2.);
end
for i=1:n1
	q(i,1) = 1.;
	for j=2:n1
		q(i,j) = xxcoll(i)^(2*j-2);
	end
end

% form c and d matrices
for j=1:n1
	ca = 2.*j-2;
	da = (2.*j-2)*(2*j+aaa-4);
	for i=1:n1
		c(i,j) = ca*xxcoll(i)^(2*j-3);
		d(i,j) = da*xxcoll(i)^(2*j-4);
	end
end

% calculate the inverse of q
qinvcoll = inv(q);

% form the aa, bb, and ww matrices
for i=1:n1
	wwcoll(i) = 0.0;
	for j=1:n1
		aacoll(i,j) = 0.;
		bbcoll(i,j) = 0.;
		for k=1:n1
			aacoll(i,j) = aacoll(i,j)+c(i,k)*qinvcoll(k,j);
			bbcoll(i,j) = bbcoll(i,j)+d(i,k)*qinvcoll(k,j);
		end
		wwcoll(i) = wwcoll(i) + z(j)*qinvcoll(j,i);
	end
end

clear c d

if kon == 0
fprintf('A matrix (aacoll) \n')
disp(aacoll)
fprintf('B matrix (bbcoll) \n')
disp(bbcoll) 
fprintf('Qinv matrix (qinvcoll) \n')
disp(qinvcoll) 
fprintf('Collocation points (xxcoll) \n')
disp(xxcoll)
fprintf('W matrix (wwcoll) \n')
disp(wwcoll)
end

		
			
% filename coll_interp.m
function y = coll_interp(ncol,dd,xb)
% ncol = number of internal collocation points, symmetric polynomials
% dd = the coefficients in the polynomial (vector)
% xb is a vector of the x values to evaluate the polynomial at

nn = size(xb);
for j=1:nn(2)
	x = xb(j);
	sum = 0;
	for i=1:ncol+1
		sum = sum*x*x + dd(ncol+2-i);
	end
	y(j) = sum;
end



function planar(nx,kon)                                                
%   this subroutine computes the matrices for collocation without       
%   symmetry, using collocation points in table 4-3.                    
%   input variable                                                     
%     nx = number of collocation points, including two end points                     
%                                                                       
%   output variables, must be declared global                                                    
%     aoc = matrix for first derivative, eq.  4-103                       
%     boc = matrix for second derivative, eq. 4-103                       
%     qinv = matrix for q inverse, eq. 4-101                               
%     xoc = vector of collocations points, from table 4-3                 
%     woc = vector of weights, eq. 4-106                                  
% 
global xoc woc aoc boc qinv                                                                      
if(nx<3) 
  nx=3;
  fprintf('\n ***************************************** \n')
  fprintf('The number of collocation points must be between \n')
  fprintf('3 and 7 (inclusive).  nx has been changed to %g\n\n',nx)
end
if(nx>7)
  nx=7;
  fprintf('\n ***************************************** \n')
  fprintf('The number of collocation points must be between \n')
  fprintf('3 and 7 (inclusive).  nx has been changed to %g\n\n',nx)
end
% now 3<=nx<=7
if kon==0 
	fprintf('The total number of collocation points is %g\n\n',nx)
end

xoc(1)=0.;
if (nx==3)
  xoc(2)=0.5;
  elseif (nx==4)
  xoc(2)=0.21132486540519;
  xoc(3)=1-xoc(2);
  elseif (nx==5)
  xoc(2)=0.11270166537926;
  xoc(3)=0.5;
  xoc(4)=1-xoc(2);
  elseif (nx==6)
  xoc(2)=0.06943184420297;
  xoc(3)=0.33000947820757;
  xoc(4)=1-xoc(3);
  xoc(5)=1-xoc(2);
  else
  xoc(2)=0.04691007703067;
  xoc(3)=0.23076534494716;
  xoc(4)=0.5;
  xoc(5)=1-xoc(3);
  xoc(6)=1-xoc(2);
end
xoc(nx)=1.;
format long                                                         
for i=1:nx,                                                        
    r(i,i) = 0.0;                                                     
    aoc(i,i) = 0.0;                                                     
    s(i) = 1.0;                                                       
    boc(i,i) = 0.0 ;                                                    
     for j=1:nx,                                                     
        if (i~=j)                                         
           r(i,j) = 1.0/(xoc(i)-xoc(j)) ;                                     
           s(i) = s(i)*r(i,j);
        end                                            
     end                                                         
     for j=1:nx,                                                     
         jx = nx-j+1;                                                   
         if (jx<j), break, end                                         
         if (jx==j) 
            aoc(i,i) = aoc(i,i)+r(i,j);
         end                           
         if (jx>j) 
            aoc(i,i) = aoc(i,i)+r(i,j)+r(i,jx);
         end                   
     end                                                        
end                                                            
for i=1:nx                                                        
   for j=1:nx                                                     
       if (i~=j)                                           
          aoc(i,j) = s(j)*r(i,j)/s(i);                                    
          boc(i,j) = 2.0*aoc(i,j)*(aoc(i,i)-r(i,j));                           
          boc(i,i) = boc(i,i)+r(i,j)*(aoc(i,i)-r(i,j)); 
       end                       
   end                                                        
end                                                           
for i=1:nx                                                        
   qinv(1,i) = s(i);                                                    
   k = 1;                                                            
   woc(i) = 0.0 ;                                                      
   for j=1:nx                                                     
      if (j~=i)                                           
         l = k;                                                         
         k = k+1;                                                       
         qinv(k,i) = qinv(l,i);                                               
         while (l>1)                                          
            m = l-1;                                                       
            qinv(l,i) = qinv(m,i)-xoc(j)*qinv(l,i);                                  
            l = m;                                                        
         end                                                                                                                  
           qinv(1,i) = -xoc(j)*qinv(1,i);
      end                                       
   end                                                         
   for j=1:nx                                                     
      woc(i) = woc(i)+qinv(j,i)/j;
   end                                     
end

if kon == 0
fprintf('A matrix (aoc) \n')
disp(aoc)
fprintf('B matrix (boc) \n')
disp(boc) 
fprintf('Q matrix (qinv) \n')
disp(qinv) 
fprintf('Collocation points (xoc) \n')
disp(xoc)
fprintf('W matrix (woc) \n')
disp(woc)
end



% filename planar_interp.m
function y = planar_interp(nx,dd,x,x0,x1)
% nx = number of collocation points, including end points
% dd = the coefficients in the polynomial (vector)
% x is a vector of the x values to evaluate the polynomial at
% The points start at x0 and go to x1, and the xb below has to 
% be scaled to go from 0 to 1.

xb = (x-x0)/(x1-x0);
nn = size(xb);
for j = 1:nn(2)
	x = xb(j);
	sum = 0;
	for i=1:nx
		sum = sum*x + dd(nx-i+1);
	end
	y(j) = sum;
end