% OC_unsteady.m
% Solves a transient diffusion problem using orthogonal collocation.
% 	need planar.m

global xoc woc aoc boc qinv nx

nx = 7
planar(nx,1)
for i=2:nx-1
	c0(i-1) = 0.;
end

tbreak = [0 .001 .01 .1 1];
for kkk=1:4
	tspan = [tbreak(kkk) tbreak(kkk+1)];
	[t cc] = ode45('OC_rhs',tspan,c0);
	nn=size(t);
	c(kkk,:) = cc(nn(1),:);
	c0(:) = c(kkk,:);
end

% add in the first and last point
for kkk=1:4
	d(kkk,1) = 1;
	for i=2:nx-1
		d(kkk,i)=c(kkk,i-1);
	end
	d(kkk,nx) = 0;
end

plot(xoc,d(4,:),'g-*',xoc,d(3,:),'m-*',xoc,d(2,:),'b-*',xoc,d(1,:),'r-*')
xlabel('x')
ylabel('c')
legend('t = 1.0','t = 0.1','t = 0.01','t = 0.001')
title('Diffusion Problem Solved using Orthogonal Collocation')
pause

x = 0:.02:1
for kkk=1:4
	pp = spline(xoc,d(kkk,:));
	e(kkk,:) = ppval(pp,x);
end

plot(xoc,d(4,:),'g*',xoc,d(3,:),'m*',xoc,d(2,:),'b*',xoc,d(1,:),'r*')
xlabel('x')
ylabel('c')
legend('t = 1.0','t = 0.1','t = 0.01','t = 0.001')
title('Diffusion Problem Solved using Spline Interpolation of Orthogonal Collocation Solution')
hold on
plot(x,e(4,:),'g-',x,e(3,:),'m-',x,e(2,:),'b-',x,e(1,:),'r-')


% OC_rhs.m
function ydot=OC_rhs(time,y)
global xoc woc aoc boc qinv nx

% add in the boundary conditions
c(1) = 1;
c(nx) = 0;

% transfer the y to the c
for i=1:nx-2
	c(i+1) = y(i);
end

% calculate the right-hand side
for i=2:nx-1
	sum = 0;
	for j = 1:nx
		sum = sum + boc(i,j)*c(j);
	end
	ydot(i-1) = sum;
end
ydot = ydot';



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