请给出一个计算系统可观测性矩阵的代码，用c++语言

以下是一个计算系统可观测性矩阵的C++代码，假设系统状态矢量为x，观测向量为y，系统动态方程为dx/dt = f(x,u)，观测方程为y=h(x)：

#include <iostream>
#include <Eigen/Dense>  //需要Eigen库支持

using namespace Eigen;
using namespace std;

int main()
{
    //定义系统状态向量维度
    const int n = 4;
    //定义观测向量维度
    const int m = 2;
    
    //定义系统状态矢量和观测矢量
    VectorXd x(n);
    VectorXd y(m);
    
    //定义系统动态方程和观测方程
    auto f = [](const VectorXd& x, double u){ return x + u * VectorXd::Ones(n); };
    auto h = [](const VectorXd& x){ return VectorXd({x(0), x(2)}); };
    
    //构造系统可观测性矩阵
    MatrixXd H(m, n);
    for(int i=0; i<m; i++)
        H.row(i) = h(x).transpose().partialDerivative(x, i);
    MatrixXd Obsv = H;
    for(int i=1; i<n; i++)
        Obsv = (H * f(x, 1.0)).colwise().stack(Obsv);
    
    //输出系统可观测性矩阵
    cout << "Observable matrix:\n" << Obsv << endl;
    
    return 0;
}

这个代码使用了Eigen库中的自动求导功能，可以自动计算观测方程对系统状态向量的偏导数，从而构造出系统可观测性矩阵。