function f = bvp_rhs(N,Nq,h);

  % Evaluate the right hand side integrals. 

  weights = gauss_weights(Nq,h);
  x0 = gauss_points(Nq,h,0);
  phi = basis_eval(Nq,h,x0); 
  
  for j=1:N-1
    xjm1 = (j-1)*h;
    xj = j*h;
    xq = gauss_points(Nq,h,xjm1);
    f1 = (pi^2)*sin((pi/3)*xq);
    f(j,1) = weights*(phi.*f1);
  end