（今天挺忙的先水水吧，代码直接放在这里了,另外附上少量的图解，以后如果有时间我会专门法一篇文章讲有限元）

下面这个是matlab平面四节点等参元的刚度矩阵计算程序

第一个是test.mlx，这个部分是计算刚度矩阵的符号矩阵的，首先对于四节点单元，其编号是逆时针顺序编号的，即：

对于四节点单元，其单元刚度矩阵是8x8的矩阵，其符号矩阵计算如下

test.mlx (这个用于测试， mlx文件只要在左边文件栏点击新建live script就行)

% 四节点单元刚度矩阵函数 E杨氏模量, nu 泊松比, h板厚度
syms x 
syms y 
syms E nu h

M = [
    [1, x, y, x*y,0,0,0,0]
    [0,0,0,0, 1, x, y, x*y]
    ];

A = [
    [1 -1 -1 1 0 0 0 0]
    [0 0 0 0 1 -1 -1 1]
    [1 1 -1 -1 0 0 0 0]
    [0 0 0 0 1 1 -1 -1]
    [1 1 1 1 0 0 0 0]
    [0 0 0 0 1 1 1 1]
    [1 -1 1 -1 0 0 0 0]
    [0 0 0 0 1 -1 1 -1]
    ];

N = M/A;  % M * inv(A) , 为形状函数矩阵

% 对x, y等参变量求导并得到对应的应变矩阵

% 对x, y等参变量求导并得到对应的应变矩阵
N_x = diff(N,x);
N_y = diff(N,y);

B = [N_x(1,:)
     N_y(2,:)
     N_y(1,:) + N_x(2,:)];

clear N_x N_y

D = E /(1-nu^2) * [
[1 nu 0]
[nu 1 0]
[0 0 (1- nu)/2]
];

k_e = transpose(B) * D * B;
Ke = h*int(int(k_e,x,[-1, 1]), y, [-1,1]);  % 积分求解对应的四节点刚度矩阵
Ke2 = Ke/(E * h/(1-nu^2));
simplify(Ke2)

求解结果是这样的(也可以使用一个.m文件来运行)

这个的正确性可以参考 《结构有限元分析》 刘永寿 第二版 第51页, 上面的刚度矩阵将a/b = 1代入之后和这个刚度矩阵结果相同。
因此我们求解单元刚度矩阵只需要将值代入就行了

下面附上可以直接复制粘贴使用的MATLAB函数代码
这个函数的调用方法是K_e = elem401(E, nu, h)
其中E是杨氏模量，nu是泊松比，h是平面板的厚度
比如elem401(2e5, 0.3, 10)其参数是杨氏模量200Gpa(由于单位是mm,而Pa是$N/m^2$ ，所以$2\times 10^{11}Pa=2\times 10^{5}N/m{m}^2)$
上面的调用方法得到刚度单位的是$N/mm$ (你使用mm为单位计算)
当然使用m为单位计算可以elem401(2e11, 0.3, 0.01)直接返回对应的刚度矩阵

elem401.m

% 四节点单元刚度矩阵函数 E杨氏模量, nu 泊松比, h板厚度
function Ke = elem401(E,nu,h)
syms x 
syms y 

M = [
    [1, x, y, x*y,0,0,0,0]
    [0,0,0,0, 1, x, y, x*y]
    ];

A = [
    [1 -1 -1 1 0 0 0 0]
    [0 0 0 0 1 -1 -1 1]
    [1 1 -1 -1 0 0 0 0]
    [0 0 0 0 1 1 -1 -1]
    [1 1 1 1 0 0 0 0]
    [0 0 0 0 1 1 1 1]
    [1 -1 1 -1 0 0 0 0]
    [0 0 0 0 1 -1 1 -1]
    ];

N = M/A;  % M * inv(A) , 为形状函数矩阵

% 对x, y等参变量求导并得到对应的应变矩阵
N_x = diff(N,x);
N_y = diff(N,y);

B = [N_x(1,:)
     N_y(2,:)
     N_y(1,:) + N_x(2,:)];

clear N_x N_y

D = E * h/(1-nu^2) * [
[1 nu 0]
[nu 1 0]
[0 0 (1- nu)/2]
];

k_e = transpose(B) * D * B;
Ke = int(int(k_e,x,[-1, 1]), y, [-1,1]);  % 积分求解对应的四节点刚度矩阵
end

另外需要注意的是，单元刚度矩阵一定是奇异矩阵，即$det(K_e)=0$， 如果需要显示矩阵，使用vpa(K_e)即可，但是注意这个会损失精度导致求解$det$不为0

下面是八节点单元的matlab代码
调用方法是K_e = Elem801(E,nu,h)
节点编号顺序为

源程序Elem801.m

function Ke = Elem801(E, nu,h)
    syms x y
    % 对于平面八节点单元，共有
    M = [
        [1 x y x^2 x*y y^2 x^2*y x*y^2 0 0 0 0 0 0 0 0]
        [0 0 0 0 0 0 0 0 1 x y x^2 x*y y^2 x^2*y x*y^2]
        ];
    
    
    % 定义八节点单元的x_e和y_e
    x_e = [-1 0 1 1 1 0 -1 -1];
    y_e = [-1 -1 -1 0 1 1 1 0];
    
    % 创建初始矩阵A
    A = zeros(16,16);
    for i = 1:8
        A(2*i-1:1:2 *i,:) = subs(M,[x,y],[x_e(i),y_e(i)]);
    end
    
    % 计算形状函数矩阵 
    N = M/A;
    
    % 使用求导方法求解对应的位移函数
    N_x = diff(N,x);
    N_y = diff(N,y);
    B = [N_x(1,:)
         N_y(2,:)
         N_y(1,:) + N_x(2,:)
         ];
    clear N_x N_y
    % 刚度矩阵D
    D = E/(1-nu^2) * [
        [1 nu 0]
        [nu 1 0]
        [0 0 (1- nu)/2]
        ];
    
    ke_pre = transpose(B)*D*B;  % 注意不能使用B'
    Ke = h * int(int(ke_pre,x,[-1, 1]), y, [-1,1]); % 积分得到刚度矩阵
    % 注意不能使用vpa函数, vpa函数会产生截断误差
    clear ke_pre
end
% Test Code: Elem801(2e5,0.3, 10)