function [ u, I ] = PDE_solver ( x, y, mu, sigma, PDE_rhs )

%*****************************************************************************80
%
%% PDE_solver() ...
%
%  Input:
%
%    x: equispaced discretization nodes.
%
%    y: random variables (vector of size 1 x N_rv).
%
%    mu: mean of the random field to be constructed.
%
%    sigma=[sigma_1 ... sigma_N_rv]: variance of each part of the random field.
%
%    PDE_rhs: function handle of the rhs of the problem, PDE_rhs = @(x) ... 
%    where x is a row vector.
%
%  Output:
%
%    u: PDE solution at x.
%
%    I: spatial integral of u, the quantity of interest.
%
  N = length(x)-1; 
  h = (x(end)-x(1))/N;

  x_mid = x(1:end-1)+h/2; 
%
%  Construct the random field.
%
  N_rv = length ( y ); 
  what_reg = @(x) max(1,ceil(x*N_rv));
    
  RF = zeros(N_rv,1);
  for i = 1:N_rv
    RF(i) = mu + sigma(i)*y(i);
  end 
  a = RF(what_reg(x_mid));
%    
%  Approximate the integrals by mid-point rule.
%
  f_mid = PDE_rhs(x_mid);
  f_int = f_mid(2:end-1);  
  rhs = h/2 * [f_int(1); f_int(1:end-1)'+f_int(2:end)'; f_int(end) ];

  A = spdiags( ...
    [ -1/h*[a(2:end-1);0], 1/h*(a(1:end-1)+a(2:end)), -1/h*[0;a(2:end-1)] ], ...
    [-1 0 1], N-1, N-1 );
%
%  Add boundary terms as first and last entries.
%
  u = [ 0; A\rhs; 0 ];
%
%  Compute the integral of u using the trapezoidal rule.
%
  I = [h/2 h*ones(1,N-1) h/2] * u; 

  return
end
