function [] = radpendel()
% Konstanten
a  = 1.00;
m1 = 0.10;
m2 = 0.15;

% Anfangs und Endzeitpunkt
t0 =   0.00;
dt =   0.01;
t1 =  10.00;

% Stationärer Punkt
x_stat = a*m1/sqrt(m2*m2-m1*m1);

% Startwerte
x0 = x_stat;
v0 = 0.1;

[T   ,Y   ] = ode45(@myFun,   [t0:dt:t1],[x0 v0]);
[Tlin,Ylin] = ode45(@myFunLin,[t0:dt:t1],[x0 v0]);

figure; plot(T, [Y(:,1) Ylin(:,1)]); 
grid on;

end



function dY = myFun(t,Y)
% Konstanten
    a  = 1.00;
    r  = 0.20;
    m1 = 0.10;
    m2 = 0.15;
    J  = 0.5*m1*r*r;
    g  = 9.81;

% Zusammengefasste Konstanten/Variablen
    m3 = m1 + J/(r*r);
    zaehler = Y(1)*Y(1)+a*a;

% Nichtlineare DGL
    dY = zeros(2,1);
    dY(1) = Y(2);
    dY(2) = m1*g...
            - m2*g*Y(1)/sqrt(zaehler)...
            - m2*a*a*Y(1)*Y(2)*Y(2)/(zaehler*zaehler);
    dY(2) = dY(2,1)/(m3 + m2*Y(1)*Y(1)/zaehler);
end



function dY = myFunLin(t,Y)
% Konstanten
    a  = 1.00;
    r  = 0.20;
    m1 = 0.10;
    m2 = 0.15;
    J  = 0.5*m1*r*r;
    g  = 9.81;

% Stationärer Punkt
    x_stat = a*m1/sqrt(m2*m2-m1*m1);

% Zusammengefasste Konstanten/Variablen
    m3 = m1 + J/(r*r);

% Linearisierte DGL
    dY = zeros(2,1);
    dY(1) = Y(2);
    dY(2) = -(g/a)*sqrt(m2*m2-m1*m1)*(1-(m1*m1)/(m2*m2))*(Y(1)-x_stat);
    dY(2) = dY(2,1)/(m3 + m1*m1/(m2));
end


