# Computing a double Integral using two Riemann sums & graphing multiple isosurfaces

I'm trying to compute the potential and electric field of a uniformly charged spherical shell and plot the results in the space outside the shell using MATLAB. I've tried everything from for loops to sum commands and more and I just can't even get the right expression for the double integral (because I can calculate that my result gives solutions that don't match analytical solutions).

From my research, it seems to best way to plot this would be to plot isosurface plots for different isovalues on the same figure to show how the potential and electric field change with distance from the sphere. But that is my second problem. I can't even compute the double integral properly. So, my question is twofold:

1. how can I compute a double integral using Riemann sums that contain symbolic variables acting as constants? I am purposely avoiding using a double integral command such as integral2() or a Reimann command like trapz.

2. how can I plot the resulting expression in 3 dimensions considering I would have a function of three variables?

The current state of my code to calculate the potential (which I define first and then later will take the gradient to get the electric field) is as follows:

syms r phi theta % symbolic values representing a random point in space
outside a uniformly charged spherical shell

% constants
N = 3; % total number of increments to divide the riemann sum into
Q = +20 .* 1e-9; % total charge on sphere
R = 1; % radius of sphere is 1 meter
sigma = Q/(4*pi*(R^2)); % uniform surface charge denisty
eps0 = 8.854e-12;
kC = 1/(4*pi*eps0); % coulomb's constant
dphi = 2*pi/N; % discretizing the interval over which I sum
dtheta = pi/N;

% Performing two riemann sums
% dA = R^2.*cos(thetaprime).*dphi.*dtheta; small area on surface of sphere
% in spherical coordinates
phiprime = linspace(0,2*pi,N);
thetaprime = linspace(0,pi,N);
dInt1 = sym(zeros(0,N)); % preallocating space allowing the presence of
symbols
dInt2 = sym(zeros(0,N));

for e = 1: length(phiprime)
for m = 1: length(thetaprime)
dInt1(m) = dInt1 +
((R^2).*cos(thetaprime).*dphi.*dtheta)./sqrt((r.*sin(theta).*cos(phi) -
R.*sin(thetaprime).*cos(phiprime)).^2 + (r.*sin(theta).*sin(phi) -
R.*sin(thetaprime).*sin(phiprime)).^2 + (r.*cos(theta) -
R.*cos(phiprime)).^2)
end

dInt2(e) = dInt2 + sum(dInt1);

end
Int = sum(Int2);

Vtot = kC*sigma.*Int; % total symbolic expression for electric potential

%% Computing Values for Vtot at various points in space and plotting them

rinterval = linspace(1,4,12);
phiinterval = linspace (0,2*pi,12);
thetainterval = linspace(0,pi,12);
[X,Y,Z] = meshgrid(rinterval,phiinterval,thetainterval);

Vcont = subs(Vtot, [r,phi,theta], {X,Y,Z}); % substituting 2D coordinates
now for symbolic values


For this code I am getting the error:

Error using symengine
Array sizes must match.

Error in sym/privBinaryOp (line 1022)