diff --git a/gpc/gpc_sample.m b/gpc/gpc_sample.m
index 17815d4c..f3f8588c 100644
--- a/gpc/gpc_sample.m
+++ b/gpc/gpc_sample.m
@@ -4,11 +4,13 @@
 %   is the number of random variables. V contains a specification of the
 %   random variables and a multiindex set (see Example). 
 %
-% Example (<a href="matlab:run_example gpc_sample">run</a>)
+% Example 1 (<a href="matlab:run_example gpc_sample 1">run</a>)
 %    I = multiindex(3,2);
 %    V = {'H', I};
 %    xi = gpc_sample(V, 5)
-%    V = {'PHp', I};
+%
+% Example 2 (<a href="matlab:run_example gpc_sample 2">run</a>)
+%    V = {'PLp', [0, 0, 0]};
 %    xi = gpc_sample(V, 3000);
 %    plot(xi(1,:), xi(2,:), '.')
 %    xlabel('Uniform'); ylabel('Gauss');
diff --git a/testing/corput.m b/testing/corput.m
deleted file mode 100644
index 86ef84e4..00000000
--- a/testing/corput.m
+++ /dev/null
@@ -1,29 +0,0 @@
-function corput
-
-clf
-xlim([0,1]);
-ylim([0,1]);
-axis equal
-x=[];
-M=100;
-for i=0:1000
-    N=(M*i):(M*(i+1)-1);
-    
-    x=[x [gb(N,3);gb(N,7)]];
-    plot(x(1,:),x(2,:),'.');
-    drawnow;
-end
-1;
-
-
-
-function g=gb( n, b )
-
-g=0;
-x=1/b;
-while n>0
-    d=mod(n,b);
-    g=g+d*x;
-    n=round((n-d)/b);
-    x=x/b;
-end
diff --git a/testing/multiindex_support.m b/testing/gpc/multiindex_support.m
similarity index 100%
rename from testing/multiindex_support.m
rename to testing/gpc/multiindex_support.m
diff --git a/testing/test_cc.m b/testing/test_cc.m
new file mode 100644
index 00000000..6d98e58b
--- /dev/null
+++ b/testing/test_cc.m
@@ -0,0 +1,22 @@
+function test_cc
+
+[x,w]=clencurt(9)
+sum(w)
+
+function [x,w]=clencurt(n)
+% Computes the Clenshaw Curtis nodes and weights
+% Adapted a code by G. von Winckel
+% http://www.scientificpython.net/1/post/2012/04/clenshaw-curtis-quadrature.html
+if n == 1
+    x = 0;
+    w = 2;
+else
+    C = zeros(n,2);
+    k = 2*(1:floor((n-1)/2));
+    C(1:2:end,1) = 2./[1, 1-k.*k];
+    C(2,2) = -(n-1);
+    V = [C; flipud(C(2:n-1,:))];
+    F = real(ifft(V)); %, n=None, axis=0))
+    x = F(1:n,2);
+    w = [F(1,1), 2*F(2:n-1,1)', F(n,1)];
+end
diff --git a/testing/test_fejer.m b/testing/test_fejer.m
new file mode 100644
index 00000000..15fc19e1
--- /dev/null
+++ b/testing/test_fejer.m
@@ -0,0 +1,52 @@
+function test_fejer
+
+
+cc_fejer(5, 1)
+cc_fejer(5, 2)
+cc_fejer(5, 3)
+sum(cc_fejer(5, 1))
+sum(cc_fejer(5, 2))
+sum(cc_fejer(5, 3))
+
+[w1,w2,wc]=fejer(16); 
+wc(end+1)=wc(1);
+sum(wc)
+length(wc)
+
+
+function w = cc_fejer(n, mode)
+[wf1,wf2,wcc] = fejer(n);
+switch mode
+    case 1
+        w = wf1';
+    case 2
+        w = wf2(2:end)';
+    case 3
+        w = [wcc; wcc(1)]';
+    otherwise
+        error('foo');
+end
+
+
+
+
+function [wf1,wf2,wcc] = fejer(n)
+%BIT Numerical Mathematics
+%March 2006, Volume 46, Issue 1, pp 195-202
+%Fast Construction of the Fejér and Clenshaw–Curtis Quadrature Rules 
+% doi: 10.1007/s10543-006-0045-4
+
+
+
+% Weights of the Fejer2, Clenshaw-Curtis and Fejer1 quadratures
+% by DFTs. Nodes: x_k = cos(k*pi/n), n>1
+N=[1:2:n-1]'; l=length(N); m=n-l; K=[0:m-1]';
+% Fejer2 nodes: k=0,1,...,n; weights: wf2, wf2_n=wf2_0=0
+v0=[2./N./(N-2); 1/N(end); zeros(m,1)];
+v2=-v0(1:end-1)-v0(end:-1:2); wf2=ifft(v2);
+%Clenshaw-Curtis nodes: k=0,1,...,n; weights: wcc, wcc_n=wcc_0
+g0=-ones(n,1); g0(1+l)=g0(1+l)+n; g0(1+m)=g0(1+m)+n;
+g=g0/(n^2-1+mod(n,2)); wcc=ifft(v2+g);
+% Fejer1 nodes: k=1/2,3/2,...,n-1/2; vector of weights: wf1
+v0=[2*exp(i*pi*K/n)./(1-4*K.^2); zeros(l+1,1)];
+v1=v0(1:end-1)+conj(v0(end:-1:2)); wf1=ifft(v1);
diff --git a/testing/multiindex2.m b/testing/wip/multiindex2.m
similarity index 100%
rename from testing/multiindex2.m
rename to testing/wip/multiindex2.m