

set(0, 'defaultaxesfontsize',24,'defaultaxeslinewidth',2.0,...
       'defaultlinelinewidth',2.0,'defaultpatchlinewidth',2.0,...
       'defaulttextfontsize',24);


%% Random Matrix Example
fprintf('\n\n Random Matrix. Hit return to continue ....\n');
pause
n = 50;
A = randn(n,n);
b = randn(n,1);
A = A'*A;        % make sure A is symmetric

[x, res, iter, flag] = mycg( A, b, [] );

figure(1)
fprintf(' CG return after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res,'k-'); 



%% 1D Laplace equation with Dirichlet BC.
% apply CG to compute an approximate solution of the FD discretization of the BVP
%
%     - y''(x) = f(x)  x in [0,1],
%       y(0) = 0,  y(1) = 0,

fprintf('\n\n 1D BVP with Dirichlet BC. Hit return to continue ....\n');
pause

% set problem data
f     = inline( '8*ones(size(x))', 'x');
yex   = inline( '4*x.*(1-x)', 'x');


% set mesh points
n   = 49;
dx  = 1/(n+1);
x   = (0:dx:1)';

% system matrix
c   = -ones(n,1);  % subdiagonal
d   = 2*ones(n,1); % diagona
e   = -ones(n,1);  % superdiagonal
A   = spdiags([c d e], -1:1, n, n);
    
% right hand side
b   = dx*dx*f(x(2:n+1));

% solve system using CG
[y, res, iter, flag] = mycg( A, b, [] );
fprintf(' CG return after %d iterations with flag = %d\n', iter,flag)

figure(1)
semilogy((0:iter),res/res(1),'k-'); 
xlabel('Iteration')
ylabel('CG Residual (Normalized)')

% plot solution
figure(2)
xx = (0:1/100:1);
plot(xx,yex(xx),'k-'); hold on
plot(x,[0;y;0],'r--'); hold off
xlabel('x')
ylabel('y(x)')
legend( 'exact BVP solution','CG solution','Location','South')
axis([-0.1,1.1,-0.1,1.1])


%% 1D Laplace equation with Neumann BC.
% apply CG to compute an approximate solution of the FD discretization of the BVP
%
%     - y''(x) = f(x)  x in [0,1],
%       y'(0) = 0,  y'(1) = 0,

fprintf('\n\n 1D BVP with Neumann BC. Hit return to continue ....\n');
pause


% set problem data
f     = inline( '4*pi^2*cos(2*pi*x)', 'x');
yex   = inline( 'cos(2*pi*x)', 'x');

f     = inline( 'cos(0.5*pi*x)', 'x'); % \int f(x) \not= 0; no solution
%yex   = inline( 'cos(0.5*pi*x)', 'x'); % PDE has no solution yex
                                        % doesn't make sense here

% set mesh points
n   = 39;
dx  = 1/(n+1);
x   = (0:dx:1)';

% generate system matrix
c   = -ones(n+2,1);                         % subdiagonal
d   = 2*ones(n+2,1); d(1) = 1; d(n+2) = 1;  % diagonal
e   = -ones(n+2,1);                         % superdiagonal
A   = spdiags([c d e], -1:1, n+2, n+2);

% right hand side
b      = dx*dx*f(x);
b(1)   = 0.5*b(1);
b(n+2) = 0.5*b(n+2);

sum(b)
y0 = 0.3*ones(n+2,1);

[y, res, iter, flag] = mycg( A, b, [], y0 );
fprintf(' CG return after %d iterations with flag = %d\n', iter,flag)

figure(1)
semilogy((0:iter),res/res(1),'k-'); 
xlabel('Iteration')
ylabel('CG Residual (Normalized)')
%print('-depsc','cg_bvp.eps')

% plot solution
figure(2)
xx = (0:1/100:1);
plot(xx,yex(xx),'k-'); hold on
plot(x,y,'r--'); hold off
xlabel('x')
ylabel('y(x)')
legend( 'exact BVP solution','CG solution','Location','North')
axis([-0.1,1.1,-1.1,1.3])
%print('-depsc','cg_laplace1d_neumann.eps')

