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

以下是一份参考代码：

% 定义计算区域
Lx = 10; % 区域长度
Ly = 10; % 区域宽度
Nx = 50; % x方向网格数
Ny = 50; % y方向网格数
dx = Lx/Nx; % x方向网格间距
dy = Ly/Ny; % y方向网格间距

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

% 定义受力情况
P = 100; % 集中荷载
fx = zeros(Nx,Ny); % x方向受力
fy = zeros(Nx,Ny); % y方向受力
fy(Nx/2,Ny/2) = -P; % 在中心施加竖向集中荷载

% 定义边界条件
u = zeros(Nx+1,Ny+1); % x方向位移
v = zeros(Nx+1,Ny+1); % y方向位移

% 构建系数矩阵和载荷向量
K = zeros((Nx+1)*(Ny+1),(Nx+1)*(Ny+1));
F = zeros((Nx+1)*(Ny+1),1);
for i=1:Nx+1
    for j=1:Ny+1
        n = (j-1)*(Nx+1)+i; % 当前节点编号
        if i==1 || i==Nx+1 || j==1 || j==Ny+1 % 边界节点
            K(n,n) = 1;
            F(n) = 0;
        else % 内部节点
            K11 = E/(1-v^2)*dy/dx; % K矩阵元素
            K22 = E/(1-v^2)*dx/dy;
            K12 = E/(1-v^2)*v*dy/dx;
            K21 = E/(1-v^2)*v*dx/dy;
            K(n,n) = K11+K22;
            K(n,n+1) = -K12;
            K(n,n-1) = -K12;
            K(n,n+Nx+1) = -K21;
            K(n,n-(Nx+1)) = -K21;
            F(n) = fx(i,j)*dx*dy + fy(i,j)*dx*dy;
        end
    end
end

% 求解位移
u_v = K\F;
for i=1:Nx+1
    for j=1:Ny+1
        n = (j-1)*(Nx+1)+i; % 当前节点编号
        u(i,j) = u_v(n);
        v(i,j) = u_v(n+(Nx+1)*(Ny+1));
    end
end

% 计算应力
sigma_x = zeros(Nx,Ny);
sigma_y = zeros(Nx,Ny);
sigma_xy = zeros(Nx,Ny);
for i=1:Nx
    for j=1:Ny
        sigma_x(i,j) = E/(1-v^2)*((v*dy/dx)*(v*(u(i,j+1)-u(i,j))-u(i+1,j)+u(i,j))+(u(i+1,j)-u(i,j))/dx);
        sigma_y(i,j) = E/(1-v^2)*((u(i,j+1)-u(i,j))/dy+(v*dx/dy)*(v*(u(i+1,j)-u(i,j))-u(i,j+1)+u(i,j)));
        sigma_xy(i,j) = E/(1+v)/2/dx*(v*(u(i+1,j+1)-u(i+1,j)-u(i,j+1)+u(i,j))+u(i,j)-u(i+1,j+1)+u(i+1,j)+u(i,j+1));
    end
end

% 绘制网格分布图
[X,Y] = meshgrid(0:dx:Lx,0:dy:Ly);
figure(1);
plot(X',Y','k',X,Y','k');
axis equal;
title('网格分布图');

% 绘制受力情况
figure(2);
quiver(X,Y,fx,fy);
axis equal;
title('受力情况');

% 绘制附加应力分布
figure(3);
contourf(X(1:end-1,1:end-1),Y(1:end-1,1:end-1),sigma_x',20,'LineStyle','none');
colorbar;
axis equal;
title('竖向应力分布');

% 绘制应力等值线
figure(4);
sigma_vm = sqrt(sigma_x.^2+sigma_y.^2-sigma_x.*sigma_y+3*sigma_xy.^2);
contour(X(1:end-1,1:end-1),Y(1:end-1,1:end-1),sigma_vm',20);
colorbar;
axis equal;
title('应力等值线');

% 绘制应力云图
figure(5);
surf(X(1:end-1,1:end-1),Y(1:end-1,1:end-1),zeros(Nx,Ny),sigma_vm','LineStyle','none');
colorbar;
axis equal;
title('应力云图');

这份代码首先定义了计算区域的大小和网格数，以及材料参数和受力情况。然后根据网格划分和边界条件构建系数矩阵和载荷向量，求解位移。接着根据位移计算应力，并绘制网格分布图、受力情况、附加应力分布、应力等值线和应力云图。