clear;
clc;
clearglobal;

//Drot = 0.01;
Drot = 1.0;
m2   = 1.0;
l    = 1.0;
h    = l/2;
r    = 0.01;
J    = 0.25*m2*r*r + (1/12)*m2*l*l;

m1   = 1.0;

FA   = 0.0;
g    = 9.81;

function f = rechteSeite(t,y)
    x1 = y(7,1);  //Weg Schlitten
   vx1 = y(8,1);
    
    x2  = y(1,1); //Weg Pendel
    y2  = y(2,1);
    phi = y(3,1);

   vx2  = y(4,1); 
   vy2  = y(5,1);
 omega  = y(6,1);

    //Zwangskräfte Fx und Fy bestimmen:
    //QQ*F=q => F = inv(QQ)*q

    q = [-FA/m1-(1/J)*h*cos(phi)*Drot*omega-h*omega*omega*sin(phi);
               -(1/J)*h*sin(phi)*Drot*omega+h*omega*omega*cos(phi)-g];

    QQ = [(1/m1)+(1/J)*h*h*cos(phi)*cos(phi)+(1/m2), (1/J)*h*h*sin(phi)*cos(phi);
          (1/J)*h*h*sin(phi)*cos(phi)              , (1/J)*h*h*sin(phi)*sin(phi) + (1/m2)];

    F = inv(QQ)*q;
    Fx = F(1);
    Fy = F(2);

    //Pendel:
    f(1,1)=x2;   //Newton x Substitutuion   x2
    f(2,1)=y2;   //Newton y Substitutuion   y2
    f(3,1)=omega;   //Euler Substitutuion      phi


    f(4,1)=-Fx/m2;   //Newton dx2/dt
    f(5,1)=-Fy/m2-g;   //Newton dy2/dt
    f(6,1)=(-Drot*omega - h*sin(phi)*Fy -h*cos(phi)*Fx)/J;   //Euler  omega

    //Schlitten:
    f(7,1)=vx1;   //Newton x Substitutuion   x1

    f(8,1)=FA/m1 + Fx/m1;   //Newton dx1/dt
endfunction

t = 0:0.01:5;
t0=0.0;
y0 = [0,0,%pi*0.5, 0,0,0,  0, 0]';

y    = ode(y0,t0,t,rechteSeite);

plot(t,y(3,:)); //phi

//plot(t,y(8,:)); //vx1
//plot(t,y(7,:)); //x1

Code 0-1: Modell für ein invertierendes Pendel (Schlitten mit Körperpendel), Datei einachser.sce.