你是一个精通计算固体力学、有限元方法底层编程以及 MATLAB 脚本开发的高级专家。
我现在需要在 MATLAB 单机环境中编写一个二维梁单元有限元分析程序，用于求解一个经典的大变形几何非线性（GNL）基准测试问题：NAFEMS Straight Cantilever GNL Benchmark。

由于这是一个包含极度大旋转和极大大挠度（屈曲及后屈曲）的悬臂梁问题，你不能使用简单的线性梁单元，必须使用能够处理任意大旋转的二维共旋梁单元（Corotational 2D Beam Element）或更新的拉格朗日格式（Updated Lagrangian Formulation），并编写增量-迭代的非线性求解器（如 Newton-Raphson 方法，结合载荷步进）。

请注意理论推导的自洽性，在此程序中，切线刚度矩阵（Tangent Stiffness Matrix）的严格定义为：内部节点力对节点位移的偏导数（KT = df_int / dd）。请确保你的内力和刚度矩阵推导遵循这一力学通用标准。

第1部分：理论检索辅助
由于大变形共旋坐标系（Corotational Formulation）的切线刚度矩阵推导极为繁琐且容易出错，在动手编写核心列式之前，请务必先使用你的联网搜索工具。
请检索并参考关于 Corotational 2D Beam Element formulation、NAFEMS Straight Cantilever GNL benchmark 以及相关大旋转有限元列式的文献、讲义或开源实现。
请基于搜索到的可靠力学理论来指导你的编程，确保单元内力和切线刚度矩阵的绝对准确与自洽。

第2部分：算例物理与几何参数
几何尺寸：直悬臂梁，长度 L = 3.2 m。横截面为正方形，边长 b = 0.1 m，高度 h = 0.1 m。
材料属性：线弹性材料，杨氏模量 E = 2.1 x 10^11 N/m^2，泊松比 v = 0（采用平面应力假设以匹配梁理论）。
边界条件：左端（x = 0）完全固定（u=0, v=0, theta=0）。
载荷条件：右端（x = 3.2）承受集中载荷。
载荷由归一化载荷因子 NCL 参数化控制，NCL 范围从 0 到 1。
总压缩轴向载荷：Fx = -3.844 x 10^6 x NCL (N)
总横向扰动载荷：Fy = -3.844 x 10^3 x NCL (N)（横向载荷较小，用于扰乱对称性以触发屈曲）。

第3部分：分析步骤要求
程序需要通过 Batch 调用的方式自动执行以下两个分析，并输出对应的图表：

分析 A：线性屈曲分析 (Linear Buckling Analysis)
在仅施加单位轴向压缩载荷（右端 Fx = -1 N）的情况下，计算线弹性刚度矩阵 KE 和几何刚度矩阵 KG。
求解特征值问题 (KE + lambda KG)Phi = 0，提取第一阶临界屈曲载荷和屈曲模态。
参考值校验：理论欧拉临界载荷约为 Pcr = 4.22 x 10^5 N。

分析 B：几何非线性大变形分析 (Geometric Non-linear Analysis)
载荷因子 NCL 从 0 逐步增加到 1（建议使用较小的载荷步，例如 dNCL = 0.01 或自适应步长）。
在每个载荷步内使用 Newton-Raphson 迭代寻根，直到残差收敛。
由于 NCL 接近 0.1 时发生屈曲，且最终状态下梁会发生极度弯曲（梁末端甚至会卷回到固定端后方，x 坐标变为负值），共旋算法必须能够处理局部坐标系超过 90 度甚至 180 度的刚体旋转。

第4部分：可视化与图表输出要求（必须通过 MATLAB 绘图复现）
程序运行结束后，必须绘制出以下三张图表（请在代码中详细配置 figure 的展示）：

图 1：梁的一阶线性屈曲模态
绘制未变形结构（直线）和一阶屈曲模态的变形形状（稍作缩放以便观察）。
在图题或图例中打印出计算得到的临界屈曲载荷因子。

图 2：末端位移 vs. 归一化载荷曲线 (Load-Displacement Curve)
横坐标：归一化载荷因子 NCL（从 0 到 1）。
纵坐标：梁右端节点的水平位移 u 和垂直位移 v（单位：m）。
曲线特征参考：在 NCL 接近 0.1 附近，曲线应出现明显的非线性转折（屈曲现象）。最终当 NCL = 1 时，水平位移 u 接近 -5.05 m，垂直位移 v 接近 -1.35 m。请将 u 和 v 画在同一张图上，并带有图例。

图 3：大变形过程及最终变形状态
以 1:1 的真实比例（axis equal）绘制大变形梁。
建议在同一张图上用不同颜色的浅色线条绘制若干个中间载荷步的变形姿态，并用深色粗线绘制 NCL = 1 时的最终变形姿态。
能够直观看到梁由于巨大的压缩载荷发生了类似“卷曲”的形态。

第5部分：自动化执行与自修复调试 (Agentic Workflow)
在给出完整的 MATLAB 源码后，请务必利用你环境中的代码执行工具（如 MATLAB 引擎接口或终端执行节点）直接运行该程序。
如果遇到任何语法报错、矩阵维度不匹配、或是 Newton-Raphson 迭代不收敛等 Bug，请仔细阅读错误日志，自行分析力学推导或代码逻辑上的原因。
请自主修改代码并重新运行，不断循环此“执行-调试-修复”过程。直到程序能够稳定计算到 NCL=1，且输出的最终位移结果（u 接近 -5.05，v 接近 -1.35）与基准参考值完全吻合，最后再将最终验证通过的正确代码呈现给我。