2D plot from slicing a surfplot

Question

I'm trying to plot a 2D slice from a surf plot I have. My code to generate the data first...

c11 = 216e9;
c12=145e9;
c44=129e9;
rho = 7900;

C=[c11,c12,c12,0,0,0;
    c12,c11,c12,0,0,0;
    c12,c12,c11,0,0,0;
    0,0,0,c44,0,0;
    0,0,0,0,c44,0;
    0,0,0,0,0,c44]; 
  
phase_vel1=zeros(181,361);
phase_vel2=zeros(181,361);
phase_vel3=zeros(181,361);
group_pvel=zeros(181,361);
group_svvel=zeros(181,361);
group_shvel=zeros(181,361);

% voight notation matrix
Voightmat=[1,5,4;5,2,6;4,6,3];

v=0;
for theta=-pi:pi/180:pi
    v=v+1;
    u=0;
    for phi=0:pi/180:pi
       u=u+1;
    n=[cos(theta)*sin(phi),sin(theta)*sin(phi),cos(phi)];
    L_11=(C(1,1)*n(1)^2+C(6,6)*n(2)^2+C(5,5)*n(3)^2+2*C(1,6)*n(1)*n(2)+2*C(1,5)*n(1)*n(3)+2*C(5,6)*n(2)*n(3));
    L_12=(C(1,6)*n(1)^2+C(2,6)*n(2)^2+C(4,5)*n(3)^2+(C(1,2)+C(6,6))*n(1)*n(2)+(C(1,4)+C(5,6))*n(1)*n(3)+(C(4,6)+C(2,5))*n(2)*n(3));
    L_13=(C(1,5)*n(1)^2+C(4,6)*n(2)^2+C(3,5)*n(3)^2+(C(1,4)+C(5,6))*n(1)*n(2)+(C(1,3)+C(5,5))*n(1)*n(3)+(C(3,6)+C(4,5))*n(2)*n(3));
    L_22=(C(6,6)*n(1)^2+C(2,2)*n(2)^2+C(4,4)*n(3)^2+2*C(2,6)*n(1)*n(2)+2*C(4,6)*n(1)*n(3)+2*C(2,4)*n(2)*n(3));
    L_23=(C(5,6)*n(1)^2+C(2,4)*n(2)^2+C(3,4)*n(3)^2+(C(4,6)+C(2,5))*n(1)*n(2)+(C(3,6)+C(4,5))*n(1)*n(3)+(C(2,3)+C(4,4))*n(2)*n(3));
    L_33=(C(5,5)*n(1)^2+C(4,4)*n(2)^2+C(3,3)*n(3)^2+2*C(4,5)*n(1)*n(2)+2*C(3,5)*n(1)*n(3)+2*C(3,4)*n(2)*n(3));
   
    Christoffel_mat=[L_11, L_12, L_13;
                    L_12,L_22,L_23;
                    L_13,L_23,L_33];   
                
    [ev,d]=eig(Christoffel_mat);
    
    % eigvecs = polarisation vectors
    pl=ev(:,3);
    psv=ev(:,1);
    psh=ev(:,2);
    
    % eigvals ~ slowness (phase)
    phase_vel1(u,v)=sqrt(rho./d(1,1)); % quasi-shear vertical
    phase_vel2(u,v)=sqrt(rho./d(2,2)); % quasi-shear horizontal
    phase_vel3(u,v)=sqrt(rho./d(3,3)); % quasi-longitudinal
    
    % calculate group vel
    GVL=zeros(1,3);
    GVSV=zeros(1,3);
    GVSH=zeros(1,3);
    
    % einstein summation
    for j=1:3
        for i=1:3
            for k=1:3
                for l=1:3
                ind1=Voightmat(i,j);
                ind2=Voightmat(k,l);
                GVL(j)=GVL(j)+C(ind1,ind2)*n(k)*pl(i)*pl(l)*phase_vel3(u,v)/rho;
                GVSV(j)=GVSV(j)+C(ind1,ind2)*n(k)*psv(i)*psv(l)*phase_vel1(u,v)/rho;
                GVSH(j)=GVSH(j)+C(ind1,ind2)*n(k)*psh(i)*psh(l)*phase_vel2(u,v)/rho;

                end
            end
        end
    end
  
    % magnitude of GV = group vel

    group_pvel(u,v)=norm(GVL);
    group_shvel(u,v)=norm(GVSH);
    group_svvel(u,v)=norm(GVSV);
    end
end

%%%%% plotting %%%%%

theta=linspace(-pi,pi,361);
phi=linspace(-pi/2,pi/2,181);
[theta1,phi1]=meshgrid(theta,phi);

[x1,y1,z1]=sph2cart(theta1,phi1,v2);
surf(x1,y1,z1,v2,'EdgeColor','none')
xlabel('x')
ylabel('y')
zlabel('z')
title('SH group')
colorbar

where x1,y1,z1 have dimension 181x361 double. I want to take the following slice

hold on
patch([-5000,-5000,5000,5000],[-4500,-4500,4500,4500],[5000,-5000,-5000,5000],'w','FaceAlpha',0.7);


% Plot an arbitrary origin axis for the slicing plane
plot3([-5000,-5000],[-4500,-4500],[-5000,5000],'r','linewidth',3);

which produces the plot

Next I want to make a meshgrid for the data to interpolate over

% Create x,y,z coordinates of the data
% z is a 181x361 double from data
[x,y]=meshgrid(1:361,1:361);
% Create x and y over the slicing plane
xq=linspace(-5000,5000,361);
yq=linspace(-4500,4500,361);

% Interpolate over the surface
zq=interp2(x,y,z,xq,yq); 

But I run into the problem when interpolating

Error using griddedInterpolant
GridVectors must define a grid whose size is compatible with the Values array.

Error in interp2>makegriddedinterp (line 226)
    F = griddedInterpolant(varargin{:});

Error in interp2 (line 134)
        F = makegriddedinterp(X, Y, V, method,extrap);

Error in spherical_slowness_curve_v2_tant_original (line 190)
zq=interp2(x,y,z,xq,yq);

How can I define the meshgrid correctly? I want to plot a 2d section of this slice...

For example, here is a plot of one slice in the y direction

%Polar plot at x deg
rho = v2(120,:);
theta = linspace(0,2*pi,length(rho));
polarplot(theta,rho)
title('Shear Group Velocity')