function [bewegungsgleichung] = build_ode
clc;
syms m1 m2 m3 l1 l2 l3 l0 Mtorq g;
syms varalpha varbeta vargamma dvaralpha dvarbeta dvargamma ddvaralpha dddvaralpha;
%Ableiteiteoperatoren
q = varalpha;
qp = dvaralpha;
qpp = ddvaralpha;
abl = [q; qp];
ablp = [qp; qpp];
%Ortsvektoren
R1 = [1/2*cos(varalpha)*l1;1/2*sin(varalpha)*l1];
R2 = [cos(varalpha)*l1+1/2*cos(varbeta)*l2;sin(varalpha)*l1+1/2*sin(varbeta)*l2];
R3 = [cos(varalpha)*l1+cos(varbeta)*l2-1/2*cos(vargamma)*l3 ;sin(varalpha)*l1+sin(varbeta)*l2-1/2*sin(vargamma)*l3];
%Geschwindigkeiten
V1 = jacobian(R1,abl)*ablp;
V2 = jacobian(R2,abl)*ablp;
V3 = jacobian(R3,abl)*ablp;
%Kinetische Energie
T = 1/2*m1*(V1.')*V1+1/24*m1*l1^2*(dvaralpha)^2+1/24*m2*l2^2*(dvarbeta)^2+1/2*m3*(V3.')*V3+1/24*m3*l3^2*(dvargamma)^2;
%Kinetische Energie: Lagrange Gleichung
T2 = jacobian(T, q);
T1 = jacobian(T, qp);
T1 = jacobian(T1, abl)*ablp;
Tlagrangestarr= T1 - T2;
%Nicht Konservative Kraft
Q = Mtorq;
%Potentielle Energie
Vpot = m1*g*[0 1]*R1 +m2*g*[0 1]*R2 +m3*g*[0 1]*R3;
%Potentielle Energie: Lagrange Gleichung
Vpot = jacobian(Vpot, q);
%Lagrange Gleichung
Lagrange = Tlagrangestarr + Vpot - Q;
%Zerlegung der Bewegungsgleichung in Ihre Einzelteile
%Masse
Mass = subs(Lagrange, {ddvaralpha, g, Mtorq}, {1, 0, 0});
%Generallisierte Kraefte
Qe = -subs(Lagrange, {ddvaralpha}, {0});
%Ueberpruefe ob Term fehlt
%--------------------Bewegungsgleichung---------------------------
bewegungsgleichung = Mass\Qe;
%dgl = matlabFunction(dx,'vars',{varalpha,l1,m1,m2,m3,g,Mtorq});