The first issue is that one of your inputs to griddata is symbolic. It must be numeric. Use the double function to convert Qs1 to doubles.
P = griddata(Lp_x11, Vs, double(Qs1), X, Z);
Next, your scatter3 function call does not use variables that exist. Variables are case sensitive. Use X, Z and P. In addition, the inputs to scatter3 must be vectors, not arrays. Use a colon to turn your arrays into vectors.
scatter3(X(:), Z(:), P(:), '.')
Surface plots (meshes, contours) will not work with the data as you currently have it since the data is more of a line than a surface. You have gridded the data, but most of the values are NaN. Scatter3 may be the best choice.
You could also try just plotting in 3D.
plot3(Lp_x11, Vs, double(Qs1))
Upon closer inspection, you'll see that you are ony solving you equation along the diagonal (Vs(i),Lp_x11(i)). This is why your results of all NaN. You are essentially trying to determine an entire surface from a single line. Instead, you want to user your loop to pick one value (say Vs(i)), and solve for all values of Lp_x11. If you capture the results correctly, Qs1 and Qs2 will be num x num arrays.
At that point, you want to use interp2, not gridded data, to increase the resolution of your data.
Here's your code modified.
num=50;
Lpmax=700e-6;
Lpmin=20e-6;
Lr=2.77e-3;
Cr=2.5e-9;
f=1/(2*pi*sqrt(Lr*Cr*0.5));
w=2*pi*f;
T=1/f;
Tdead=T*0.011;
Lp_x11=linspace(Lpmin,Lpmax,num);
Vs=linspace(140,340,num);
Vo=180;
Re_x11=200;
theta=0;
% solve the D
syms D11
theta_dch=(2*pi*Tdead)/T;
for r=1:length(Vs)
X=vpasolve(sin(D11*pi)/(2*(1-D11))-Vo/Vs(r));
D_ch=X;
A2=(1./(Re_x11.*(pi.^2).*(1-D_ch))).*2.*w.*Lp_x11.*sin(pi.*D_ch).*sin(pi.*(1-D_ch));
Qlr1=(2.*Vs(r).*sin(pi.*D_ch)./(Re_x11.*pi.*(1-D_ch).*w)).*(sin(-pi.*D_ch-theta_dch)-sin(-pi.*D_ch));
Qlr2=(2.*Vs(r).*sin(pi.*D_ch)./(Re_x11.*pi.*(1-D_ch).*w)).*(sin(pi*D_ch-theta_dch)-sin(pi.*D_ch));
Qlp_S1=-(1./w).*(Vs(r).*theta_dch./(w.*Lp_x11)).*(A2./(1-D_ch)-pi.*D_ch+theta_dch./2);
Qlp_S2=(1./w).*(D_ch.*Vs(r).*theta_dch)./(w.*Lp_x11.*(1-D_ch)).*(A2./D_ch+pi.*(1-D_ch)-theta_dch./2);
% Vs will be the y axis, which corresponds to the rows of the matrix
% Capture the entire row of data using the indexing below.
Qs1(r,:)=Qlp_S1+Qlr1;
Qs2(r,:)=Qlr2+Qlp_S2;
end
xv = linspace(min(Lp_x11), max(Lp_x11), 1000);
zv = linspace (min (Vs), max (Vs), 1000);
[X,Z] = meshgrid(xv, zv);
% P = griddata(Lp_x11, Vs, double(Qs1), X, Z);
P = interp2(Lp_x11, Vs, double(Qs1), X, Z);
figure
scatter3(X(:), Z(:), P(:), '.')
grid on
figure
meshc(X, Z, P)
grid on
figure
contourf(X, Z, P)
grid on
