matrix df(double t,matrix Y,matrix ud1,matrix ud2)
{
    double c= , r=,m,r0,g;
    double s=M_PI*r*r;
    double Dx=0.5*c*r0*s*abs(Y(1))*(Y(1));
    double Dy= 0.5 *c*r0*s*abs(Y(3))*(Y(3));
    double Fmx=M_PI*r0*Y(3)*m2d(ud2)*pow(r,3);
    double Fmy=M_PI*r0*Y(1)*m2d(ud2)*pow(r,3);

    matrix dy;
    dy(0)=Y(1);
    dy(1)=(-dy-Fmx)/m;
    dy(2)=Y(3);
    dy(3)=(-m*g-dy-Fmy)/m;
}