function y=Ggen(krnl,t,knots,order)
%Ggen  Generate the matrix
%
%             G(i,j) = int_knots(1)^knots(n) krnl(t(i),s) B_j(s) ds
%
%      for i = 1, ... , length(t), j = 1, ... , n,
%      where n = length(knots) - order
%      krnl is a string variable containing the name of the kernel
%           function.  Usage of krnl is  krnl(s),  where the first
%           argument is the global variable TOUT
%
global GKERNEL TOUT SPDATA
GKERNEL = krnl;
lt = length(t);
n = length(knots) - order;
y = zeros(lt,n);
coef = 1;
for j = 1:n
   fprintf( 'Working on column (basis function) %g\n', j )
   SPDATA = spmak( knots(j:j+order), coef );
   a = knots(j);
   b = knots(j+order);
   for i = 1:lt
      TOUT = t(i);
      y(i,j) = quad('integrand',a,b);
   end
   fprintf( '\n' )
end