How to compute a double Integral using two Riemann sums and graph f(x,y,z) = c

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 doesn't give solutions that match analytical solutions).

Basically, I'm trying to compute the potential and field of a uniformly charged spherical shell and plot the results in the space outside the shell. From my research, it seems to best way to plot this would be to plot isosurface plots for different isovalues (graph f(x,y,z) = c for multiple values of c) 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?

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

Right now I am currently getting the error:

Error using symengine
Array sizes must match.

Error in sym/privBinaryOp (line 1022)
            Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

Error in  +  (line 7)
X = privBinaryOp(A, B, 'symobj::zipWithImplicitExpansion', '_plus');

Error in Numerical_Integration_of_2D_Surfaces (line 27)
      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)

Thanks in advance!