%Vorlaufwerte für ABM-Verfahren mit klassischem Runge-Kutta-V. berechnen

    %benötigt:
    %b...Balkenanzahl
    %n...Anzahl der Zeitschritte
    %h...Schrittweite
    %M_sys...Systemmassenmatrix
    %K_sys...Systemmsteifigkeitsmatrix
    %Gamma_sys...Systemdämpfungsmatrix
    %F_sys...Systemlastvektor

    %Startwerte       
    x01=zeros(3*(b+1),1);   %= x(t=0); [mm]
    x02=zeros(3*(b+1),1);   %= x°(t=0); [mm/s]
    %x02(2:3:end,1)=2;
    
    k1=zeros(3*(b+1),n+1,4);
    k2=zeros(3*(b+1),n+1,4);
    dk1=zeros(3*(b+1),n+1,4);
    dk2=zeros(3*(b+1),n+1,4);
    k1(:,1,1)=x01;
    k2(:,1,1)=x02;
    
    for i=1:3       %i...zu berechnende Zeitschritte
        %1. Hilfswert (:,:,1) -> k1
        dk1(:,i,1)=k2(:,i,1);
        dk2(:,i,1)=M_sys\(-Gamma_sys*k2(:,i,1)-K_sys*k1(:,i,1)+F_sys(:,i));
        %Funktionswert zu k1
        k1(:,i,2)=k1(:,i,1)+(h/2)*dk1(:,i,1);
        k2(:,i,2)=k2(:,i,1)+(h/2)*dk2(:,i,1);
        
        %2. Hilfswert (:,:,2) -> k2
        dk1(:,i,2)=k2(:,i,2);
        dk2(:,i,2)=M_sys\(-Gamma_sys*k2(:,i,2)-K_sys*k1(:,i,2)+F_sys(:,i));
        %Fkt.wert zu k2
        k1(:,i,3)=k1(:,i,1)+(h/2)*dk1(:,i,2);
        k2(:,i,3)=k2(:,i,1)+(h/2)*dk2(:,i,2);
        
        %3. Hilfswert (:,:,3) -> k3
        dk1(:,i,3)=k2(:,i,3);
        dk2(:,i,3)=M_sys\(-Gamma_sys*k2(:,i,3)-K_sys*k1(:,i,3)+F_sys(:,i));
        %Fkt.wert zu k3
        k1(:,i,4)=k1(:,i,1)+h*dk1(:,i,3);
        k2(:,i,4)=k2(:,i,1)+h*dk2(:,i,3);
        
        %4. Hilfswert (:,:,4) -> k4
        dk1(:,i,4)=k2(:,i,4);
        dk2(:,i,4)=M_sys\(-Gamma_sys*k2(:,i,4)-K_sys*k1(:,i,4)+F_sys(:,i));
        %Fkt.wert zu k4
        k1(:,i+1,1)=k1(:,i,1)+(h/6)*(dk1(:,i,1)+2*dk1(:,i,2)+2*dk1(:,i,3)+dk1(:,i,4));
        k2(:,i+1,1)=k2(:,i,1)+(h/6)*(dk2(:,i,1)+2*dk2(:,i,2)+2*dk2(:,i,3)+dk2(:,i,4));
    end
    k1(:,:,1);
    k2(:,:,1);