beta  = 1.0;
a     = 0.1;

rlo = 1.0;
rhi = 5.0;
nr  = 4;
dr  = (rhi - rlo) / nr;

thetalo = 0;
thetahi = pi;
ntheta  = 4;
dtheta  = (thetahi - thetalo) / ntheta;

[rs, thetas] = meshgrid(rlo:dr:rhi, thetalo:dtheta:thetahi);
[nr, ntheta] = size(rs);
Xs = zeros(nr, ntheta);
Is = zeros(nr, ntheta);
for j = 1 : ntheta
	for i = 1 : nr
		r     = rs(i, j);
		theta = thetas(i, j);
		X     = v2exact_scalar(r,    theta, beta, a);
		I     = v2i           (r, r, theta, beta, a);
		Xs(i, j) = X;
		Is(i, j) = I;

		string = sprintf(...
	'r = %5.2f theta = %6.2f int = %10.7f exact = %10.7f relerr = %9.2e', ...
			r, theta/pi*180, I, X, (I-X)/X);
		disp(string);
	end
	disp(' ');
end
