写一段matlab代码，用于计算竖向集中荷载作用在土体顶部中央时土中的附加应力，并且顶部为自由位移，其他三边位移为0。划分好网格之后画出沿深度方向网格分布图以及受力情况。要点：不需要调用代码包，而且以弹性的方式计算，将其作为平面应变问题并且包含边界的处理，计算区域的选择，并且画出计算区域以及受力情况，画出沿深度方向附加应力分布，应力等值线，应力云图

好的，我可以帮你编写这段MATLAB代码。基于你的要求，我们可以将问题转化为平面应变问题，并采用弹性理论来计算土体中的附加应力。下面是一份可能的代码：

% 定义输入参数
L = 5; % 计算区域的长度
W = 5; % 计算区域的宽度
H = 2; % 土体的高度
E = 10e6; % 土体的弹性模量
v = 0.3; % 土体的泊松比
q = 1000; % 集中荷载大小
Nx = 50; % x方向网格数
Ny = 50; % y方向网格数
Nz = 10; % z方向网格数

% 定义计算区域
x = linspace(0, L, Nx);
y = linspace(0, W, Ny);
z = linspace(0, H, Nz);
[X, Y, Z] = meshgrid(x, y, z);

% 计算附加应力
nu = 1 - 2*v;
k = E/(3*(1-2*v));
r = sqrt((X-L/2).^2 + (Y-W/2).^2);
theta = atan2(Y-W/2, X-L/2);
sigma_x = q/2 + k*nu./(H-Z).^2.*(3*W/2 - r.*cos(theta) - nu*(L/2 - r.*sin(theta)));
sigma_y = q/2 + k*nu./(H-Z).^2.*(3*L/2 - r.*sin(theta) - nu*(W/2 - r.*cos(theta)));
tau_xy = -k./(H-Z).^2.*r.*sin(theta);

% 处理边界
sigma_x(:,1,:) = 0;
sigma_x(:,end,:) = 0;
sigma_y(1,:,:) = 0;
sigma_y(end,:,:) = 0;
tau_xy(:,1,:) = 0;
tau_xy(:,end,:) = 0;
tau_xy(1,:,:) = 0;
tau_xy(end,:,:) = 0;

% 画出网格分布图
figure;
imagesc(z, y, zeros(Ny, Nz));
hold on;
plot3(zeros(Ny*Nz,1), Y(:), Z(:), 'k.', 'MarkerSize', 10);
axis equal;
title('沿深度方向网格分布图');

% 画出受力情况
figure;
quiver3(X, Y, Z, sigma_x, sigma_y, zeros(size(sigma_x)), 'AutoScale', 'off');
hold on;
plot3(L/2, W/2, H, 'ro', 'MarkerSize', 10);
axis equal;
title('受力情况');

% 画出附加应力分布
figure;
imagesc(z, y, squeeze(sqrt(sigma_x(:, :, end).^2 + sigma_y(:, :, end).^2 + 2*tau_xy(:, :, end).^2))));
hold on;
plot3(zeros(Ny*Nz,1), Y(:), Z(:,end), 'k.', 'MarkerSize', 10);
colorbar;
axis equal;
title('沿深度方向附加应力分布');

% 画出应力等值线
figure;
contour(y, z, squeeze(sqrt(sigma_x(:, :, end).^2 + sigma_y(:, :, end).^2 + 2*tau_xy(:, :, end).^2))), 20);
hold on;
plot(W/2, H, 'ro', 'MarkerSize', 10);
colorbar;
axis equal;
title('应力等值线');

% 画出应力云图
figure;
quiver3(X, Y, Z, sigma_x, sigma_y, zeros(size(sigma_x)), 'AutoScale', 'off');
hold on;
contour(y, z, squeeze(sqrt(sigma_x(:, :, end).^2 + sigma_y(:, :, end).^2 + 2*tau_xy(:, :, end).^2))), 20);
plot3(L/2, W/2, H, 'ro', 'MarkerSize', 10);
colorbar;
axis equal;
title('应力云图');

这段代码中，我们首先定义了输入参数，包括计算区域的大小、土体的弹性模量、泊松比和集中荷载大小。然后，我们使用meshgrid函数创建了计算区域的网格，并计算了每个网格点处的附加应力。接着，我们处理了边界条件，并画出了沿深度方向的网格分布图和受力情况。最后，我们画出了沿深度方向的附加应力分布、应力等值线和应力云图。

希望这份代码对你有所帮助！