基于以上代码，修改并计算沿深度方向的力与位移分布，并添加绘制xz平面上的位移大小云图，应力等值线图以及剪应力分布图。

为了沿深度方向计算力与位移分布，需要添加一个z轴坐标，并在计算过程中引入z维度的参量。以下是修改后的代码：

% 定义模型参数和常量
E = 20e6; % 基岩弹性模量
v = 0.3; % 泊松比
a = 1; % 圆形荷载半径
P = 10000; % 竖向集中荷载
nx = 50; ny = 50; nz = 50; % X方向、Y方向和Z方向的网格数目

% 定义计算区域
xmin = -5; xmax = 5; ymin = -5; ymax = 5; zmin = 0; zmax = 10;
x = linspace(xmin, xmax, nx);
y = linspace(ymin, ymax, ny);
z = linspace(zmin, zmax, nz);
[X, Y, Z] = meshgrid(x, y, z);

% 计算相对位移、边界条件和附加应力
r = sqrt(X.^2 + Y.^2);
u_r = P*a^2./(4*E)*(1-v^2)./(r.^2+a^2-2*a*r.*cos(pi/2));
u_theta = P*a^2./(4*E)*2*v./(r.^2+a^2-2*a*r.*cos(pi/2));
sigma_rr = P*a^2./(2*pi*E)*(r.^2-a^2)./((r.^2+a^2-2*a*r.*cos(pi/2)).^2);
sigma_theta_theta = P*a^2./(2*pi*E)*((r.^2+a^2)./(r.^2-a^2).*...
    (1+v) - 2*v*cos(2*atan(Y./X)))./((r.^2+a^2-2*a*r.*cos(pi/2)).^2);
sigma_xy = P*a^2./(2*pi*E)*2*v.*sin(2*atan(Y./X))./...
    ((r.^2+a^2-2*a*r.*cos(pi/2)).^2);
u_x = u_r.*X./r - u_theta.*Y./r;
u_y = u_r.*Y./r + u_theta.*X./r;
u_z = P*z/(E*pi*a^2);

% 绘制网格分布图和应力等值线、云图
figure(1)
subplot(2,2,1)
plot(x,y,'k-',x',y','k-')
title('计算区域')
xlabel('x (m)'); ylabel('y (m)')

subplot(2,2,2)
contourf(X(:,:,round(nz/2)),Z(:,:,round(ny/2)),squeeze(u_x(:,:,round(ny/2))),20)
colorbar
title('xz平面位移大小云图')
xlabel('x (m)'); ylabel('z (m)')

subplot(2,2,3)
contour(X(:,:,round(nz/2)),Z(:,:,round(ny/2)),squeeze(sigma_rr(:,:,round(ny/2))),20,'k-')
hold on
contour(X(:,:,round(nz/2)),Z(:,:,round(ny/2)),squeeze(sigma_theta_theta(:,:,round(ny/2))),20,'k--')
title('xz平面应力等值线图')
xlabel('x (m)'); ylabel('z (m)');
legend('\sigma_{rr}','\sigma_{\theta\theta}')

subplot(2,2,4)
quiver(X(:,:,round(nz/2)),Z(:,:,round(ny/2)),squeeze(sigma_xy(:,:,round(ny/2))),zeros(size(squeeze(sigma_xy(:,:,round(ny/2))))))
title('xz平面剪应力分布图')
xlabel('x (m)'); ylabel('z (m)')

% 绘制位移大小云图、应力等值线图和剪应力分布图
figure(2)
subplot(2,2,1)
contourf(X(:,:,round(ny/2)),Y(:,:,round(nz/2)),squeeze(sqrt(u_x(:,:,round(ny/2)).^2+u_y(:,:,round(ny/2)).^2)),20)
colorbar
title('xy平面位移大小云图')
xlabel('x (m)'); ylabel('y (m)')

subplot(2,2,2)
contour(X(:,:,round(ny/2)),Y(:,:,round(nz/2)),squeeze(sigma_rr(:,:,round(ny/2))),20,'k-')
hold on
contour(X(:,:,round(ny/2)),Y(:,:,round(nz/2)),squeeze(sigma_theta_theta(:,:,round(ny/2))),20,'k--')
title('xy平面应力等值线图')
xlabel('x (m)'); ylabel('y (m)');
legend('\sigma_{rr}','\sigma_{\theta\theta}')

subplot(2,2,3)
quiver(X(:,:,round(ny/2)),Y(:,:,round(nz/2)),squeeze(sigma_xy(:,:,round(ny/2))),zeros(size(squeeze(sigma_xy(:,:,round(ny/2))))))
title('xy平面剪应力分布图')
xlabel('x (m)'); ylabel('y (m)')

subplot(2,2,4)
plot(z, squeeze(u_z(:,round(nx/2))))
title('沿深度方向的位移分布图')
xlabel('z (m)'); ylabel('u_z (m)')

其中，添加了一个Z轴坐标并定义了网格数目，然后计算了沿Z轴方向的位移大小和应力分布，以及沿深度方向的位移分布，并绘制了xz平面和xy平面上的位移大小云图、应力等值线图和剪应力分布图。