l=1.0;
b=0.6;
c=0.5*b;
m=10.0;
J=(1/12)*m*(l*l+b*b);
D=1;

global l,b,c,m,J,D;

function f = rechteSeite(t,y)
//Reihenfolge y: xx,vxx,yy,vyy,phi,omega
    xx    = y(1,1); //Schwerpunt in I
    vxx   = y(2,1);
    yy    = y(3,1); //Schwerpunt in I
    vyy   = y(4,1);
    phi   = y(5,1); //Drehung von Quader in I
    omega = y(6,1);

    TT = [cos(phi),-sin(phi);sin(phi),cos(phi)];
    dTT = [-omega*sin(phi),-omega*cos(phi);omega*cos(phi),-omega*sin(phi)];

    oPPstrich = [0;c];
    oQQstrich = [0;-c];

    oPP = [xx;yy] + TT*oPPstrich; //Ortsvektor nach Punkt P in I
    oQQ = [xx;yy] + TT*oQQstrich; //Ortsvektor nach Punkt Q in I

    voPP = [vxx;vyy] + dTT*oPPstrich; //Ortsvektor-Geschw. von Punkt P in I
    voQQ = [vxx;vyy] + dTT*oQQstrich; //Ortsvektor-Geschw. von Punkt Q in I

    rPP = TT*oPPstrich; //Verbindungsvektor von Punkt S nach P in I
    rQQ = TT*oQQstrich; //Verbindungsvektor von Punkt S nach Q in I


//Vorwärts in X-Richtung:
//    skalarF1 = 10;
//    skalarF2 = 10;

//Drehen:
    skalarF1 = -10;
    skalarF2 = 10;

    F1 = TT*[skalarF1;0];
    F2 = TT*[skalarF2;0];

    FR1 = -D*voPP;
    FR2 = -D*voQQ;

    M1 = rPP(1,1)*F1(2,1) - rPP(2,1)*F1(1,1); //rPP x F1
    M2 = rQQ(1,1)*F2(2,1) - rQQ(2,1)*F2(1,1); //rQQ x F2

    MR1 = rPP(1,1)*FR1(2,1) - rPP(2,1)*FR1(1,1); //rPP x FR1
    MR2 = rQQ(1,1)*FR2(2,1) - rQQ(2,1)*FR2(1,1); //rQQ x FR2

    f(1,1) =  vxx;
    f(2,1) =  ( F1(1,1) + F2(1,1) + FR1(1,1) + FR2(1,1) )/m ;
    f(3,1) =  vyy;
    f(4,1) =  ( F1(2,1) + F2(2,1) + FR1(2,1) + FR2(2,1) )/m ;
    f(5,1) =  omega;
    f(6,1) =  ( M1 + M2 + MR1 + MR2 )/J ;
endfunction

t = 0:0.01:30;
y0 = [0,0,0,0,0,0]';
t0 = 0;
y  = ode(y0,t0,t,rechteSeite);

//plot(t,y(2,:)');
//plot(t,y(1,:)');
//plot(y(1,:)',y(3,:)');

plot(t,y(6,:)');

Code 0-1: Scilab-Code