function [uh,A,b]=oned_diff_eg(N)
% C.E. Powell, MATH64082
% Function to generate realisations of the RF a(x,omega)
% in Example 1.1 and finite difference approximations to corresponding
% solution u(x, omega) to the BVP (1.1). See Figure 1. 
%
% N = no. intervals to use in finite difference grid.
% For example, run using:
% >> [uh,A,b]=oned_diff_eg(64);

P= 50; mu=1; sigma=4; 

% set-up 1d FD mesh
% Assuming zero Dirichlet bcs and f=1;
% N= no. elements/intervals in finite difference grid

h=(1/N);  xc=(h/2):h:1-(h/2); xc=xc';
% sample a(x,omega) at element centres
a_cent=mu.*ones(N,1);   % set the mean of a
xi=-1+2.*rand(P,1);     % draw sample of the xi_k
% construct the approximate realisation of a 
 for i=1:P
     a_cent=a_cent+sigma.*((i.*pi).^(-2)).*cos(pi.*i.*xc).*xi(i);
 end
 
% Assemble finite difference linear system + solve.
v1=a_cent(1:end-1)*(1/h^(2)); v3=a_cent(2:end)*(1/h^(2));
v2=(v1+v3);
A = spdiags([-v3 v2 -v1], -1:1, N-1, N-1);
b = ones(N-1,1);

% solve linear system for interior degrees of freedom;
u_int=A\b; uh=[0;u_int;0]; 
xx=0:h:1;
figure(10);
h=subplot(1,2,1); plot(xc,a_cent); hold on; axis([0,1,0,2]); axis square;
h=subplot(1,2,2); plot(xx',uh,'-'); hold on; axis square; 


