## Code for Figure 9

### main.m

clear('all'); close('all');
n = 10;
[x, A, B, Q] = colloc(n);
a = 10;

nplot = 100;
xplot = linspace(0, 1, nplot)';

y = sin(2*pi*xplot);

Y = sin(2*pi*x);

plot(xplot, y,  x, Y, 'o')

table1 = [xplot, y];
table2 = [x, Y];

save sinfcn_1.dat table1
save sinfcn_2.dat table2


### /export/home/jbraw/courses/cbe255/content/util/matlab/colloc.m

function [r, a, b, q] = colloc (n, varargin)

nargs = nargin;

if (nargs < 1 || nargs > 3)
error ('usage: [r, a, b, q] = colloc (n, ''left'', ''right'')');
end

if (~ (isnumeric (n) && numel (n) == 1 && round (n) == n))
error ('colloc: first argument must be an integer scalar');
end

if (isnan (n) || isinf (n))
error ('colloc: NaN is invalid as N');
end

if (n < 0)
error ('colloc: first argument must be non-negative');
end

n0 = 0;
n1 = 0;

for i = 2:nargs
s = varargin{i-1};
if (~ ischar (s))
error ('colloc: expecting character string argument');
end

s = lower (s);
if (strcmp (s, 'r') || strcmp (s, 'right'))
n1 = 1;
elseif (strcmp (s, 'l') || strcmp (s, 'left'))
n0 = 1;
else
error ('colloc: unrecognized argument');
end
end

nt = n + n0 + n1;

if (nt < 1)
error ('colloc: the total number of roots must be positive');
end

alpha = 0;
beta = 0;

%% Compute roots.

[dif1, dif2, dif3, r] = jcobi (n, n0, n1, alpha, beta);

%% First derivative weights.

a = zeros (nt, nt);
for i = 1:nt
a(i,:) = dfopr (n, n0, n1, i, 1, dif1, dif2, dif3, r)';
end

%% Second derivative weights.

b = zeros (nt, nt);
for i = 1:nt
b(i,:) = dfopr (n, n0, n1, i, 2, dif1, dif2, dif3, r)';
end

id = 3;
q = dfopr (n, n0, n1, 0, id, dif1, dif2, dif3, r);

end

%% The following routines (JCOBI, DIF, DFOPR, INTRP, AND RADAU)
%% are the same as found in Villadsen, J. and M.L. Michelsen,
%% Solution of Differential Equation Models by Polynomial
%% Approximation, Prentice-Hall (1978) pages 418-420.
%%
%% Cosmetic changes (elimination of arithmetic IF statements, most
%% GO TO statements, and indentation of program blocks) made by:
%%
%% John W. Eaton
%% Department of Chemical Engineering
%% The University of Texas at Austin
%% Austin, Texas 78712
%%
%% June 6, 1987
%%
%%
%% Further cosmetic changes made August 20, 1987
%%
%% Translated from Fortran December 14, 2006

function vect = dfopr (n, n0, n1, i, id, dif1, dif2, dif3, root)

%%   Villadsen and Michelsen, pages 133-134, 419
%%
%%   Input parameters:
%%
%%     N      : The degree of the Jacobi polynomial, (i.e. the number
%%              of interior interpolation points)
%%
%%     N0     : Determines whether x = 0 is included as an
%%              interpolation point
%%
%%                n0 = 0  ==>  x = 0 is not included
%%                n0 = 1  ==>  x = 0 is included
%%
%%     N1     : Determines whether x = 1 is included as an
%%              interpolation point
%%
%%                n1 = 0  ==>  x = 1 is not included
%%                n1 = 1  ==>  x = 1 is included
%%
%%     I      : The index of the node for which the weights are to be
%%              calculated
%%
%%     ID     : Indicator
%%
%%                id = 1  ==>  first derivative weights are computed
%%                id = 2  ==>  second derivative weights are computed
%%                id = 3  ==>  Gaussian weights are computed (in this
%%                             case, I is not used).
%%
%%   Output parameters:
%%
%%     DIF1   : vector containing the first derivative
%%              of the node polynomial at the zeros
%%
%%     DIF2   : vector containing the second derivative
%%              of the node polynomial at the zeros
%%
%%     DIF3   : vector containing the third derivative
%%              of the node polynomial at the zeros
%%
%%     VECT   : vector of computed weights

if (n0 ~= 0 && n0 ~= 1)
error ('dfopr: n0 not equal to 0 or 1');
end

if (n1 ~= 0 && n1 ~= 1)
error ('dfopr: n1 not equal to 0 or 1');
end

if (n < 0)
error ('dfopr: n less than 0');
end

nt = n + n0 + n1;

if (id ~= 1 && id ~= 2 && id ~= 3)
error ('dfopr: id not equal to 1, 2, or 3');
end

if (id ~= 3)
if (i < 1)
error ('dfopr: index less than zero');
end

if (i > nt)
error ('dfopr: index greater than number of interpolation points');
end
end

if (nt < 1)
error ('dfopr: number of interpolation points less than 1');
end

%% Evaluate discretization matrices and Gaussian quadrature
%% weights.  Quadrature weights are normalized to sum to one.

vect = zeros (nt, 1);

if (id ~= 3)
for j = 1:nt
if (j == i)
if (id == 1)
vect(i) = dif2(i)/dif1(i)/2.0;
else
vect(i) = dif3(i)/dif1(i)/3.0;
end
else
y = root(i) - root(j);
vect(j) = dif1(i)/dif1(j)/y;
if (id == 2)
vect(j) = vect(j)*(dif2(i)/dif1(i) - 2.0/y);
end
end
end
else
y = 0.0;

for j = 1:nt

x = root(j);
ax = x*(1.0 - x);

if (n0 == 0)
ax = ax/x/x;
end

if (n1 == 0)
ax = ax/(1.0 - x)/(1.0 - x);
end

vect(j) = ax/dif1(j)^2;
y = y + vect(j);

end

vect = vect/y;

end

end

function [dif1, dif2, dif3, root] = jcobi (n, n0, n1, alpha, beta)

%%   Villadsen and Michelsen, pages 131-132, 418
%%
%%   This subroutine computes the zeros of the Jacobi polynomial
%%
%%      (ALPHA,BETA)
%%     P  (X)
%%      N
%%
%%   Input parameters:
%%
%%     N      : The degree of the Jacobi polynomial, (i.e. the number
%%              of interior interpolation points)
%%
%%     N0     : Determines whether x = 0 is included as an
%%              interpolation point
%%
%%                N0 = 0  ==>  x = 0 is not included
%%                N0 = 1  ==>  x = 0 is included
%%
%%     N1     : Determines whether x = 1 is included as an
%%              interpolation point
%%
%%                N1 = 0  ==>  x = 1 is not included
%%                N1 = 1  ==>  x = 1 is included
%%
%%     ALPHA  : The value of alpha in the description of the jacobi
%%              polynomial
%%
%%     BETA   : The value of beta in the description of the Jacobi
%%              polynomial
%%
%%     For a more complete explanation of alpha and beta, see Villadsen
%%     and Michelsen, pages 57 to 59
%%
%%   Output parameters:
%%
%%     ROOT   : vector containing the n + n0 + n1 zeros of the node
%%              polynomial used in the interpolation routine
%%
%%     DIF1   : vector containing the first derivative
%%              of the node polynomial at the zeros
%%
%%     DIF2   : vector containing the second derivative
%%              of the node polynomial at the zeros
%%
%%     DIF3   : vector containing the third derivative
%%              of the node polynomial at the zeros

if (n0 ~= 0 && n0 ~= 1)
error ('jcobi: n0 not equal to 0 or 1');
end

if (n1 ~= 0 && n1 ~= 1)
error ('jcobi: n1 not equal to 0 or 1');
end

if (n < 0)
error ('jcobi: n less than 0');
end

nt = n + n0 + n1;

if (nt < 1)
error ('jcobi: number of interpolation points less than 1');
end

dif1 = zeros (nt, 1);
dif2 = zeros (nt, 1);
dif3 = zeros (nt, 1);
root = zeros (nt, 1);

%% First evaluation of coefficients in recursion formulas.
%% recursion coefficients are stored in dif1 and dif2.

ab = alpha + beta;
ap = beta*alpha;
dif1(1) = (ad/(ab + 2.0) + 1.0)/2.0;
dif2(1) = 0.0;

if (n >= 2)
for i = 2:n

z1 = i - 1.0;
z = ab + 2*z1;
dif1(i) = (ab*ad/z/(z + 2.0) + 1.0)/2.0;

if (i == 2)
dif2(i) = (ab + ap + z1)/z/z/(z + 1.0);
else
z = z*z;
y = z1*(ab + z1);
y = y*(ap + y);
dif2(i) = y/z/(z - 1.0);
end

end
end

%% Root determination by newton method with suppression of
%% previously determined roots.

x = 0.0;

for i = 1:n

done = false;

while (~ done)

xd = 0.0;
xn = 1.0;
xd1 = 0.0;
xn1 = 0.0;

for j = 1:n
xp = (dif1(j) - x)*xn  - dif2(j)*xd;
xp1 = (dif1(j) - x)*xn1 - dif2(j)*xd1 - xn;
xd = xn;
xd1 = xn1;
xn = xp;
xn1 = xp1;
end

zc = 1.0;
z = xn/xn1;

if (i ~= 1)
for j = 2:i
zc = zc - z/(x - root(j-1));
end
end

z = z/zc;
x = x - z;

if (abs(z) <= 1.0e-09)
done = true;
end

end

root(i) = x;
x = x + 0.0001;

end

%% Add interpolation points at x = 0 and/or x = 1.

nt = n + n0 + n1;

if (n0 ~= 0)
for i = 1:n
j = n + 1 - i;
root(j+1) = root(j);
end
root(1) = 0.0;
end

if (n1 == 1)
root(nt) = 1.0;
end

%% Use recursion formulas to evaluate derivatives of node polynomial
%%
%%                   N0     (ALPHA,BETA)           N1
%%     P  (X)  =  (X)   *  P (X)         *  (1 - X)
%%      NT                   N
%%
%% at the interpolation points.

for i = 1:nt
x = root(i);
dif1(i) = 1.0;
dif2(i) = 0.0;
dif3(i) = 0.0;
for j = 1:nt
if (j ~= i)
y = x - root(j);
dif3(i) = y*dif3(i) + 3.0*dif2(i);
dif2(i) = y*dif2(i) + 2.0*dif1(i);
dif1(i) = y*dif1(i);
end
end
end

end


## Data files

### sinfcn_1.dat

# Created by Octave 3.8.1, Thu Aug 21 10:29:15 2014 CDT
# name: table1
# type: matrix
# rows: 100
# columns: 2
0 0
0.0101010101010101 0.06342391965656451
0.0202020202020202 0.1265924535737493
0.0303030303030303 0.1892512443604102
0.04040404040404041 0.2511479871810792
0.05050505050505051 0.3120334456984871
0.06060606060606061 0.3716624556603275
0.07070707070707072 0.4297949120891717
0.08080808080808081 0.4861967361004687
0.09090909090909091 0.5406408174555976
0.101010101010101 0.5929079290546405
0.1111111111111111 0.6427876096865394
0.1212121212121212 0.6900790114821119
0.1313131313131313 0.7345917086575333
0.1414141414141414 0.7761464642917569
0.1515151515151515 0.8145759520503357
0.1616161616161616 0.8497254299495144
0.1717171717171717 0.8814533634475821
0.1818181818181818 0.9096319953545183
0.1919191919191919 0.9341478602651067
0.202020202020202 0.954902241444074
0.2121212121212121 0.9718115683235417
0.2222222222222222 0.9848077530122081
0.2323232323232323 0.9938384644612541
0.2424242424242424 0.998867339183008
0.2525252525252525 0.9998741276738751
0.2626262626262627 0.9968547759519424
0.2727272727272728 0.9898214418809327
0.2828282828282829 0.9788024462147786
0.2929292929292929 0.9638421585599422
0.303030303030303 0.9450008187146685
0.3131313131313131 0.9223542941045814
0.3232323232323233 0.8959937742913359
0.3333333333333334 0.8660254037844385
0.3434343434343435 0.8325698546347712
0.3535353535353536 0.795761840530832
0.3636363636363636 0.7557495743542583
0.3737373737373738 0.7126941713788627
0.3838383838383839 0.6667690005162917
0.393939393939394 0.6181589862206051
0.4040404040404041 0.5670598638627704
0.4141414141414142 0.5136773915734063
0.4242424242424243 0.4582265217274105
0.4343434343434344 0.4009305354066136
0.4444444444444445 0.3420201433256685
0.4545454545454546 0.2817325568414297
0.4646464646464647 0.2203105327865404
0.4747474747474748 0.1580013959733494
0.4848484848484849 0.09505604330418288
0.494949494949495 0.03172793349806766
0.5050505050505051 -0.03172793349806786
0.5151515151515152 -0.09505604330418307
0.5252525252525253 -0.1580013959733501
0.5353535353535354 -0.2203105327865406
0.5454545454545455 -0.2817325568414298
0.5555555555555556 -0.3420201433256687
0.5656565656565657 -0.4009305354066142
0.5757575757575758 -0.4582265217274107
0.5858585858585859 -0.5136773915734061
0.595959595959596 -0.567059863862771
0.6060606060606061 -0.6181589862206053
0.6161616161616162 -0.6667690005162918
0.6262626262626263 -0.7126941713788629
0.6363636363636365 -0.7557495743542587
0.6464646464646465 -0.7957618405308321
0.6565656565656566 -0.8325698546347713
0.6666666666666667 -0.8660254037844388
0.6767676767676768 -0.895993774291336
0.686868686868687 -0.9223542941045817
0.696969696969697 -0.9450008187146683
0.7070707070707072 -0.9638421585599422
0.7171717171717172 -0.9788024462147787
0.7272727272727273 -0.9898214418809327
0.7373737373737375 -0.9968547759519424
0.7474747474747475 -0.9998741276738751
0.7575757575757577 -0.998867339183008
0.7676767676767677 -0.9938384644612541
0.7777777777777778 -0.9848077530122081
0.787878787878788 -0.9718115683235417
0.797979797979798 -0.9549022414440739
0.8080808080808082 -0.9341478602651064
0.8181818181818182 -0.9096319953545182
0.8282828282828284 -0.881453363447582
0.8383838383838385 -0.8497254299495144
0.8484848484848485 -0.8145759520503358
0.8585858585858587 -0.7761464642917566
0.8686868686868687 -0.7345917086575332
0.8787878787878789 -0.6900790114821114
0.888888888888889 -0.6427876096865389
0.8989898989898991 -0.5929079290546402
0.9090909090909092 -0.5406408174555974
0.9191919191919192 -0.4861967361004688
0.9292929292929294 -0.4297949120891711
0.9393939393939394 -0.3716624556603272
0.9494949494949496 -0.3120334456984862
0.9595959595959597 -0.2511479871810794
0.9696969696969697 -0.1892512443604106
0.9797979797979799 -0.126592453573749
0.9898989898989899 -0.06342391965656452
1 -2.449293598294706e-16



### sinfcn_2.dat

# Created by Octave 3.8.1, Thu Aug 21 10:29:15 2014 CDT
# name: table2
# type: matrix
# rows: 10
# columns: 2
0.01304673574141415 0.08188327832197055
0.06746831665550772 0.4113329152832688
0.1602952158504878 0.8453203753408128
0.2833023029353764 0.9781881186306984
0.4255628305091844 0.4508367402262003
0.5744371694908156 -0.4508367402262001
0.7166976970646236 -0.9781881186306984
0.8397047841495122 -0.8453203753408127
0.9325316833444922 -0.4113329152832692
0.9869532642585859 -0.08188327832197104