How to plot streamlines of 3D vector field?

I want to plot the streamlines of the 3D vector field of the electric field of a static charge distribution. The field is calculated using Coulomb's law. The code works fine, it gives the desired vector field (using quiver3, no problem there). However, it fails to compute the streamlines of the vector field. When I run the code below, I get the error message: "Error using gridedInterpolant. Interpolation requires at least two sample points in each dimension." I don't know what I am doing wrong, the rest of the code runs fine. When I did it for the 2D field, I had to transpose my matrices so the code for 2D is streamline(x',y',Ex',Ey',x',y'). Help would be much appreciated!

function p=vecEfield3D(a,b,c,PX,PY,PZ,Q,d)
% The function p=vecEfield3D(a,b,c,PX,PY,PZ,Q,d) returns the E-field of the
% static charge distribution specified by the arrays PX,PY,PZ and Q. The
% parameters a,b,c and d allow to specify the grid.
 
 
set(0,'DefaultTextInterpreter','latex');
set(0,'DefaultAxesFontSize',15);
set(0,'DefaultTextFontSize',15);
set(0,'DefaultLegendInterpreter','latex');
 
 
[x,y,z]=ndgrid(a:d:-a,b:d:-b,c:d:-c);
X=linspace(a,-a,length(x));
Y=linspace(b,-b,length(y));
Z=linspace(c,-c,length(z));
 
K=8.9875e9;
 
RX=zeros(length(X),length(Y),length(Z),length(Q));
RY=zeros(length(X),length(Y),length(Z),length(Q));
RZ=zeros(length(X),length(Y),length(Z),length(Q));
ME=zeros(length(X),length(Y),length(Z),length(Q));
EX=zeros(length(X),length(Y),length(Z),length(Q));
EY=zeros(length(X),length(Y),length(Z),length(Q));
EZ=zeros(length(X),length(Y),length(Z),length(Q));
 
for n=1:length(Q)
    for i=1:length(X)
        for j=1:length(Y)
            for l=1:length(Z)
                RX(i,j,l,n)=X(i)-PX(n);
                RY(i,j,l,n)=Y(j)-PY(n);
                RZ(i,j,l,n)=Z(l)-PZ(n);
                if (X(i)==PX(n))&&(Y(j)==PY(n))&&(Z(l)==PZ(n))
                    ME(i,j,l,n)=0;
                else
                    ME(i,j,l,n)=K*Q(n)/(norm([RX(i,j,l,n),RY(i,j,l,n),RZ(i,j,l,n)]))^2;
                end
                
                EX(i,j,l,n)=ME(i,j,l,n)*RX(i,j,l,n)/(norm([RX(i,j,l,n),RY(i,j,l,n),RZ(i,j,l,n)]));
                EY(i,j,l,n)=ME(i,j,l,n)*RY(i,j,l,n)/(norm([RX(i,j,l,n),RY(i,j,l,n),RZ(i,j,l,n)]));
                EZ(i,j,l,n)=ME(i,j,l,n)*RZ(i,j,l,n)/(norm([RX(i,j,l,n),RY(i,j,l,n),RZ(i,j,l,n)]));
                
            end
        end
    end
end
 
Ex=sum(EX,4);
Ey=sum(EY,4);
Ez=sum(EZ,4);
 
p=figure;
set(gcf,'position',[10000 10000 10000 400]);
subplot(1,2,1);
    box on
    hold on
    quiver3(x,y,z,Ex,Ey,Ez);
    plot3(PX,PY,PZ,'r.','MarkerSize',20);
    xlabel('$x$-Direction [m]');
    ylabel('$y$-Direction [m]');
    zlabel('$z$-Direction [m]');
    title('\textbf{Electrostatic vector field}');
    view(3);
subplot(1,2,2);
    box on
    hold on
    streamline(x,y,z,Ex,Ey,Ez,x,y,z);
    plot3(PX,PY,PZ,'r.','MarkerSize',20);
    xlabel('$x$-Direction [m]');
    ylabel('$y$-Direction [m]');
    zlabel('$z$-Direction [m]');
    title('\textbf{Electrostatic vector field lines}');
    view(3);
hold off