beta  = 1.0;
a     = 0.1;

rlo = 0.3;
rhi = 5.0;
nr  = 50;
dr  = (rhi - rlo) / nr;

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

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

surf(rs, thetas/pi*180, Xs);
shading interp;
