帮我写一个计算拉格朗日方程组的matlab代码

以下是计算拉格朗日方程组的 MATLAB 代码：

syms t q1(t) q2(t) q1dot(t) q2dot(t) m1 m2 g l1 l2

% 定义拉格朗日函数
L = 1/2*m1*(l1*q1dot)^2 + 1/2*m2*((l1*q1dot)^2 + (l2*q2dot)^2 + 2*l1*l2*q1dot*q2dot*cos(q1-q2)) - m1*g*l1*cos(q1) - m2*g*(l1*cos(q1) + l2*cos(q2));

% 计算拉格朗日方程组
eq1 = diff(diff(L, q1dot), t) - diff(L, q1) == 0;
eq2 = diff(diff(L, q2dot), t) - diff(L, q2) == 0;

% 将符号变量转换成函数句柄
f1 = matlabFunction(eq1, 'vars', [t q1 q2 q1dot q2dot]);
f2 = matlabFunction(eq2, 'vars', [t q1 q2 q1dot q2dot]);

% 定义初始条件和时间范围
q1_0 = pi/4;
q2_0 = pi/4;
q1dot_0 = 0;
q2dot_0 = 0;
tspan = [0 10];

% 求解微分方程
[t, y] = ode45(@(t, y) [q1dot(t); q2dot(t); f1(t, y(1), y(2), y(3), y(4)); f2(t, y(1), y(2), y(3), y(4))], tspan, [q1_0 q2_0 q1dot_0 q2dot_0]);

% 绘制图像
plot(t, y(:,1), t, y(:,2));
legend('q1', 'q2');
xlabel('t');
ylabel('q');

这段代码使用符号计算工具箱中的 `syms` 命令定义了拉格朗日函数，然后使用 `diff` 命令计算出拉格朗日方程组的两个方程，再将符号变量转换成函数句柄，使用 MATLAB 的 `ode45` 命令求解微分方程，并最终绘制出系统的运动轨迹。