It looks like streamline does not accepts coordinates coming from your spherical construction. If you start form a cartesian grid it works (note the filtering of the field to avoid the singularities)
%VECTOR FIELD AND FIELD LINES OF A STATIC ELECTRIC DIPOLE%
%(1) Go ahead and define the constant 1/(4*pi*epsilon0)
epsilon0 = 8.85*((10)^(-12));
constant1 = ((4*pi)^(-1))*epsilon0;
%(2) Define the grid over which the vector field will be defined. We want a thickened spherical shell that "just barely"
% excludes the origin to avoid singularities. Keep in mind that in regards to phi(the azimuthal angle) 0 = 2*pi; We don't want redundant
%points that will cause the interpolation feature when drawing field lines to produce a "non-unique points" error.
%Same with theta(the elevation angle) which is awkward given the usual spherical coordinate convention, but for some reason only
%half the field shows up if we only vary theta from 0 to pi. By letting theta vary from 0 to 2*pi, hopefully we don't run into "non-unique
%points" error.
%(n-1) will represent the actual number of points in each vector
n=25;
rmin = 0.2;
rmax = 2.0;
%%%%% THIS IS THE NEW PART %%%%%%
x = linspace(-rmax,rmax,n);
[X,Y,Z] = meshgrid(x,x,x);
% field
Fx = X;
Fy = Y;
Fz = Z;
% remove singularities
idx = sqrt(X.^2+Y.^2+Z.^2) < rmin;
Fx(idx) = NaN;
Fy(idx) = NaN;
Fz(idx) = NaN;
%%%%%% END %%%%%%
%(5) Create a new meshgrid for the innermost sphere, this will serve as our starting points for drawing streamlines. Remember we must
%also convert this Spherical meshgrid to Cartesian.
phi2 = linspace(0,2*pi - ((n)^(-1))*2*pi,n-1);
theta2 = linspace(0,2*pi - ((n)^(-1))*2*pi,n-1);
rstart = rmin;
[Phi2,Theta2,Rstart] = meshgrid(phi2,theta2,rstart);
[Xstart,Ystart,Zstart] = sph2cart(Phi2,Theta2,Rstart);
%(7) From the documentation of streamlines it seems that startx, starty, startz must be column vectors consisting of all the x,y,z starting points
% respectively. It seems that corresponding entries must in said column vectors form starting points.
%Extract vectors from [Xstart,Ystart,Zstart] and then transpose
Xstart2 = Xstart(:);
Ystart2 = Ystart(:);
Zstart2 = Zstart(:);
Xstart3 = Xstart2';
Ystart3 = Ystart2';
Zstart3 = Zstart2';
%(6) Plot the "position vector field" from (4) and try to draw streamlines. Hopefully we do not get the error "grid vectors must contain
%unique points". If we succesfully plot the vector field and field lines then comment this section out and move on to the actual field of
%the dipole.
x = linspace(-1,1,n);
[X,Y,Z] = meshgrid(x,x,x);
figure
quiver3(X,Y,Z,Fx,Fy,Fz)
hold on
figure, axis equal
streamline(X,Y,Z,Fx,Fy,Fz,Xstart3,Ystart3,Zstart3)
view([1 1 1]);
