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
%fprintf('The total number of collocation points is %g\n\n',nx)

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
