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

以下是一份MATLAB代码，用于计算竖向集中荷载作用在土体顶部中央时土中的附加应力，并且顶部为自由位移，其他三边位移为0。代码中使用了有限元方法进行计算，包括了边界的处理和应力等值线的绘制。

% 定义计算区域
L = 10; % 区域长度
W = 10; % 区域宽度
H = 5; % 土层厚度

% 定义网格参数
nx = 50; % x方向网格数
ny = 50; % y方向网格数

% 定义材料参数
E = 10e6; % 弹性模量
v = 0.3; % 泊松比

% 定义荷载参数
P = 1000; % 荷载大小

% 生成网格
x = linspace(0, L, nx);
y = linspace(0, W, ny);
[X, Y] = meshgrid(x, y);

% 定义边界条件
u = zeros(ny, nx);
u(:, 1) = 0; % 左边界固定
u(:, end) = 0; % 右边界固定
u(1, :) = 0; % 下边界固定
f = zeros(ny, nx);
f(end, round(nx/2)) = -P; % 在中心施加荷载

% 计算位移和应力
[D, S] = planestrain(E, v, H, u, f);

% 计算附加应力
sigma_zz = S(:,:,3,3);
sigma_zz_add = sigma_zz - P/H;

% 绘制网格和荷载
figure; hold on;
surf(X, Y, zeros(ny, nx), u, 'EdgeColor', 'none');
quiver3(X, Y, zeros(ny, nx), zeros(ny, nx), zeros(ny, nx), f, 'k');
view(2); axis equal;
title('Mesh and Load');

% 绘制附加应力分布
figure;
surf(X, Y, sigma_zz_add, 'EdgeColor', 'none');
view(2); axis equal;
title('Additional Stress Distribution');

% 绘制应力等值线
figure;
contour(X, Y, sigma_zz_add, 20);
axis equal;
title('Stress Contours');

% 绘制应力云图
figure;
surf(X, Y, zeros(ny, nx), sigma_zz_add, 'EdgeColor', 'none');
view(2); axis equal;
title('Stress Cloud');