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

好的，以下是我为您编写的MATLAB代码，用于计算竖向集中荷载作用在土体顶部中央时土中的附加应力，并绘制相关图形：

% 定义计算区域大小和网格数目
L = 10; % 长度
W = 10; % 宽度
nx = 100; % x方向网格数目
ny = 100; % y方向网格数目

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

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

% 计算单元尺寸
dx = L/nx;
dy = W/ny;

% 定义节点坐标
x = linspace(0, L, nx+1);
y = linspace(0, W, ny+1);
[X, Y] = meshgrid(x, y);

% 定义节点编号
node_nums = reshape(1:numel(X), size(X));

% 定义边界条件
bc = zeros(numel(X), 1);
bc(node_nums(:,1)) = 1; % 左边界
bc(node_nums(:,end)) = 1; % 右边界
bc(node_nums(1,:)) = 1; % 上边界

% 定义刚度矩阵和荷载向量
K = zeros(numel(X));
F = zeros(numel(X), 1);

for i = 1:nx
    for j = 1:ny
        if bc(node_nums(i,j)) == 0 % 内部节点
            % 定义单元节点编号
            n1 = node_nums(i,j);
            n2 = node_nums(i+1,j);
            n3 = node_nums(i+1,j+1);
            n4 = node_nums(i,j+1);
            
            % 计算单元刚度矩阵
            D = E/(1-v^2)*[1 v 0; v 1 0; 0 0 (1-v)/2];
            Ke = (D/dx/dy)*[1 -1 1 -1; -1 1 -1 1; 1 -1 1 -1; -1 1 -1 1];
            
            % 组装整个系统的刚度矩阵和荷载向量
            K([n1 n2 n3 n4],[n1 n2 n3 n4]) = K([n1 n2 n3 n4],[n1 n2 n3 n4]) + Ke;
            F(n1) = F(n1) - P/dx/dy/2;
            F(n2) = F(n2) - P/dx/dy/2;
            F(n3) = F(n3) - P/dx/dy/2;
            F(n4) = F(n4) - P/dx/dy/2;
        end
    end
end

% 处理边界条件
K(bc==1,:) = 0;
K(:,bc==1) = 0;
K(bc==1,bc==1) = eye(sum(bc==1));
F(bc==1) = 0;

% 求解位移和应力
u = K\F;
D = E/(1-v^2)*[1 v 0; v 1 0; 0 0 (1-v)/2];
Sx = zeros(size(X));
Sy = zeros(size(Y));
Sxy = zeros(size(X));
for i = 1:nx
    for j = 1:ny
        if bc(node_nums(i,j)) == 0 % 内部节点
            % 定义单元节点编号
            n1 = node_nums(i,j);
            n2 = node_nums(i+1,j);
            n3 = node_nums(i+1,j+1);
            n4 = node_nums(i,j+1);
            
            % 计算位移梯度
            B = 1/4/dx*[y(j+1)-y(j) y(j)-y(j) y(j)-y(j+1) y(j+1)-y(j+1); x(i+1)-x(i) x(i)-x(i) x(i)-x(i+1) x(i+1)-x(i+1)];
            
            % 计算应力
            Si = D*B*u([n1 n2 n3 n4])';
            Sx(i,j) = Si(1);
            Sy(i,j) = Si(2);
            Sxy(i,j) = Si(3);
        end
    end
end

% 绘制网格分布图
figure;
hold on;
for i = 1:size(X,1)
    plot(X(i,:), Y(i,:), 'k');
end
for j = 1:size(X,2)
    plot(X(:,j), Y(:,j), 'k');
end
axis equal;
axis([0 L 0 W]);
title('网格分布图');

% 绘制受力情况
figure;
hold on;
quiver(X, Y, Sx, Sy, 1, 'LineWidth', 1);
axis equal;
axis([0 L 0 W]);
title('受力情况');

% 绘制附加应力分布
figure;
hold on;
contourf(X, Y, Sxy, 20);
colorbar;
axis equal;
axis([0 L 0 W]);
title('附加应力分布');

% 绘制应力等值线
figure;
hold on;
contour(X, Y, sqrt(Sx.^2 + Sy.^2), 20);
colorbar;
axis equal;
axis([0 L 0 W]);
title('应力等值线');

% 绘制应力云图
figure;
hold on;
surf(X, Y, sqrt(Sx.^2 + Sy.^2));
shading interp;
colorbar;
axis equal;
axis([0 L 0 W]);
title('应力云图');