% OZDEMO
% This program creates the Ornstein-Zernike example in Chapter 3.
% [H,C]=OZDEMO returns the solution on a grid with a mesh
% spacing of 1/256.
%
function [h,c]=ozdemo
global L U rho
n=257;
epsilon=.1; sigma=2.0; rho=.2; beta=10; L=9;
dx=L/(n-1); r=0:dx:L; r=r'; 
% 
% Compute the potential and store it in a global variable.
%
U=elj(r,sigma,epsilon,beta);
%
tol=[1.d-8,1.d-8];
x=zeros(2*n,1);
parms=[40,80,-.1];
[sol, it_hist, ierr] = nsoli(x,'oz',tol);
%
% Unpack h and c.
%
h=sol(1:n); c=sol(n+1:2*n);
%
% Plot the solution.
%
figure(1)
subplot(1,2,1)
plot(r,h,'-');
ylabel('h','Rotation',1);
xlabel('r');
subplot(1,2,2)
plot(r,c,'-');
ylabel('c','Rotation',1);
xlabel('r');
%
% Do a second solve with constant forcing terms.
% Are you getting the same results each time? 
%
fb=it_hist(1,1);
parms=[40,80,-.1]; x=zeros(2*n,1);
[sola, it_hist1, ierr] = nsoli(x,'oz',tol,parms);
norm(sola-sol)
%
% plot residual vs iteration counter
%
figure(2)
ni=length(it_hist(:,1));
n1=length(it_hist1(:,1));
semilogy(0:ni-1,it_hist(:,1)/fb,'-',...
0:n1-1,it_hist1(:,1)/fb,'--');
legend('default','.1');
xlabel('Nonlinear iterations');
ylabel('Relative residual');
%
% plot residual vs function counter
%
figure(3)
semilogy(it_hist(:,2),it_hist(:,1)/fb,'-',...
it_hist1(:,2),it_hist1(:,1)/fb,'--');
xlabel('Function evaluations');
ylabel('Relative residual');
legend('default','.1');

%
function u=elj(r,sigma,epsilon,beta)
n2=length(r);
ra=r(2:n2);
r12=(sigma./ra).^12; r6=(sigma./ra).^6;
ua=exp(-4*beta*epsilon*(r12-r6));
u=[0,ua']';