    % % % % % % % % % % % % % % % % % % % % % % % % % %
    %                                                 %
    % Programmieraufgabe SS07, A29, DGL               %
    %                                                 %
    % % % % % % % % % % % % % % % % % % % % % % % % % %

    % Variablenbezeichnungen:
    %
    % n : Dimension des Problems
    % m : Anzahl der Schritte 
    % h : Schrittweite
    % y : abhaengige Variablen (Dimension n x m+1), je nach Verfahren;
    %     y_euler : Loesung mit Euler Verfahren
    %     y_veuler: Loesung mit verbesstertem Euler Verfahren
    %     y_rk    : Loesung mit Runge Kutta Verfahren
    % x : unabhaengige Variablen (Ortsvariable x oder Zeit t)
    % x0: Startwert der unabhaengigen Variablen
    % xe: Endwert der unabhaengigen Variablen




    % Aufgabe (a): 
    % ------------
    close all 
    clear all



    n = 1; 
    x0 = 0; 
    xe = 1; 

    m=1;
    h = (xe-x0)/m; 
    y_euler = zeros(m,1);
    y_euler(1) = 1;
    y_veuler = zeros(m,1);
    y_veuler(1) = 1;
    y_rk = zeros(m,1);
    y_rk(1) = 1;
    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(it+1) = euler(@fa,n,x,h,y_euler(it));
        y_veuler(it+1) = veuler(@fa,n,x,h,y_veuler(it));
        y_rk(it+1) = rungekutta(@fa,n,x,h,y_rk(it));
    end

    h1=h;
    lsg1_euler=y_euler;
    lsg1_veuler=y_veuler;
    lsg1_rk=y_rk;


    m=2;
    h = (xe-x0)/m; 
    y_euler = zeros(m,1);
    y_euler(1) = 1;
    y_veuler = zeros(m,1);
    y_veuler(1) = 1;
    y_rk = zeros(m,1);
    y_rk(1) = 1;
    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(it+1) = euler(@fa,n,x,h,y_euler(it));
        y_veuler(it+1) = veuler(@fa,n,x,h,y_veuler(it));
        y_rk(it+1) = rungekutta(@fa,n,x,h,y_rk(it));
    end

    h2=h;
    lsg2_euler=y_euler;
    lsg2_veuler=y_veuler;
    lsg2_rk=y_rk;

    m=4;
    h = (xe-x0)/m; 
    y_euler = zeros(m,1);
    y_euler(1) = 1;
    y_veuler = zeros(m,1);
    y_veuler(1) = 1;
    y_rk = zeros(m,1);
    y_rk(1) = 1;
    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(it+1) = euler(@fa,n,x,h,y_euler(it));
        y_veuler(it+1) = veuler(@fa,n,x,h,y_veuler(it));
        y_rk(it+1) = rungekutta(@fa,n,x,h,y_rk(it));
    end

    h3=h;
    lsg3_euler=y_euler;
    lsg3_veuler=y_veuler;
    lsg3_rk=y_rk;


    m=8;
    h = (xe-x0)/m; 
    y_euler = zeros(m,1);
    y_euler(1) = 1;
    y_veuler = zeros(m,1);
    y_veuler(1) = 1;
    y_rk = zeros(m,1);
    y_rk(1) = 1;
    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(it+1) = euler(@fa,n,x,h,y_euler(it));
        y_veuler(it+1) = veuler(@fa,n,x,h,y_veuler(it));
        y_rk(it+1) = rungekutta(@fa,n,x,h,y_rk(it));
    end

    h4=h;
    lsg4_euler=y_euler;
    lsg4_veuler=y_veuler;
    lsg4_rk=y_rk;


    lsg_ex(1)=2*exp(x0)-(1+x0);
    for it=1:m
       x=x0+it*h;
       lsg_ex(it+1)=2*exp(x)-(1+x);
    end


    scrsz = get(0,'ScreenSize');
    figure('Position',[0 500 1400 500])
    subplot(1,3,1)
    plot(x0:h1:xe,lsg1_euler,'k')
    hold on
    plot(x0:h2:xe,lsg2_euler,'g')
    plot(x0:h3:xe,lsg3_euler,'m')
    plot(x0:h4:xe,lsg4_euler,'b')
    plot(x0:h4:xe,lsg_ex,'r')
    hold off
    xlabel('x')
    ylabel('y(x)')
    axis square
    title('AWA 1, Euler-Verfahren')
    subplot(1,3,2)
    plot(x0:h1:xe,lsg1_veuler,'k')
    hold on
    plot(x0:h2:xe,lsg2_veuler,'g')
    plot(x0:h3:xe,lsg3_veuler,'m')
    plot(x0:h4:xe,lsg4_veuler,'b')
    plot(x0:h4:xe,lsg_ex,'r')
    hold off
    xlabel('x')
    ylabel('y(x)')
    axis square
    title('AWA 1, verbessertes Euler-Verfahren')
    subplot(1,3,3)
    plot(x0:h1:xe,lsg1_rk,'k')
    hold on
    plot(x0:h2:xe,lsg2_rk,'g')
    plot(x0:h3:xe,lsg3_rk,'m')
    plot(x0:h4:xe,lsg4_rk,'b')
    plot(x0:h4:xe,lsg_ex,'r')
    hold off
    axis square
    xlabel('x')
    ylabel('y(x)')
    legend('m=1','m=2','m=4','m=8','exakt',2)
    title('AWA 1, Runge Kutta Verfahren')



    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %


    % Aufgabe (b): 
    % ------------

    n = 2; 
    x0 = 0; 
    xe = 5; 

    m=20;
    h = (xe-x0)/m; 
    y_euler = zeros(n,m+1);
    y_euler(1,1) = 3;
    y_euler(2,1) = 1;
    y_veuler = zeros(n,m+1);
    y_veuler(1,1) = 3;
    y_veuler(2,1) = 1;
    y_rk = zeros(n,m+1);
    y_rk(1,1) = 3;
    y_rk(2,1) = 1;

    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(:,it+1) = euler(@fb,n,x,h,y_euler(:,it));
        y_veuler(:,it+1) = veuler(@fb,n,x,h,y_veuler(:,it));
        y_rk(:,it+1) = rungekutta(@fb,n,x,h,y_rk(:,it));
    end

    h1=h;
    lsg1_euler=y_euler;
    lsg1_veuler=y_veuler;
    lsg1_rk=y_rk;


    m=50;
    h = (xe-x0)/m; 
    y_euler = zeros(n,m+1);
    y_euler(1,1) = 3;
    y_euler(2,1) = 1;
    y_veuler = zeros(n,m+1);
    y_veuler(1,1) = 3;
    y_veuler(2,1) = 1;
    y_rk = zeros(n,m+1);
    y_rk(1,1) = 3;
    y_rk(2,1) = 1;

    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(:,it+1) = euler(@fb,n,x,h,y_euler(:,it));
        y_veuler(:,it+1) = veuler(@fb,n,x,h,y_veuler(:,it));
        y_rk(:,it+1) = rungekutta(@fb,n,x,h,y_rk(:,it));
    end

    h2=h;
    lsg2_euler=y_euler;
    lsg2_veuler=y_veuler;
    lsg2_rk=y_rk;


    m=100;
    h = (xe-x0)/m; 
    y_euler = zeros(n,m+1);
    y_euler(1,1) = 3;
    y_euler(2,1) = 1;
    y_veuler = zeros(n,m+1);
    y_veuler(1,1) = 3;
    y_veuler(2,1) = 1;
    y_rk = zeros(n,m+1);
    y_rk(1,1) = 3;
    y_rk(2,1) = 1;

    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(:,it+1) = euler(@fb,n,x,h,y_euler(:,it));
        y_veuler(:,it+1) = veuler(@fb,n,x,h,y_veuler(:,it));
        y_rk(:,it+1) = rungekutta(@fb,n,x,h,y_rk(:,it));
    end

    h3=h;
    lsg3_euler=y_euler;
    lsg3_veuler=y_veuler;
    lsg3_rk=y_rk;

    scrsz = get(0,'ScreenSize');
    figure('Position',[0 500 1400 500])
    subplot(1,3,1)
    plot(x0:h1:xe,lsg1_euler(1,:),'k')
    hold on
    plot(x0:h1:xe,lsg1_euler(2,:),'k--')
    plot(x0:h2:xe,lsg2_euler(1,:),'g')
    plot(x0:h2:xe,lsg2_euler(2,:),'g--')
    plot(x0:h3:xe,lsg3_euler(1,:),'r')
    plot(x0:h3:xe,lsg3_euler(2,:),'r--')
    hold off
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    title('AWA 2, Raeuber-Beute-Modell, Euler-Verfahren')
    subplot(1,3,2)
    plot(x0:h1:xe,lsg1_veuler(1,:),'k')
    hold on
    plot(x0:h1:xe,lsg1_veuler(2,:),'k--')
    plot(x0:h2:xe,lsg2_veuler(1,:),'g')
    plot(x0:h2:xe,lsg2_veuler(2,:),'g--')
    plot(x0:h3:xe,lsg3_veuler(1,:),'r')
    plot(x0:h3:xe,lsg3_veuler(2,:),'r--')
    hold off
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    title('AWA 2, Raeuber-Beute-Modell, verbessertes Euler-Verfahren')
    subplot(1,3,3)
    plot(x0:h1:xe,lsg1_rk(1,:),'k')
    hold on
    plot(x0:h1:xe,lsg1_rk(2,:),'k--')
    plot(x0:h2:xe,lsg2_rk(1,:),'g')
    plot(x0:h2:xe,lsg2_rk(2,:),'g--')
    plot(x0:h3:xe,lsg3_rk(1,:),'r')
    plot(x0:h3:xe,lsg3_rk(2,:),'r--')
    hold off
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    legend('m=20, x(t)','m=20, y(t)','m=50, x(t)','m=50, y(t)','m=100, x(t)','m=100, y(t)',2)
    title('AWA 2, Raeuber-Beute-Modell, Runge Kutta Verfahren')


    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %


    % Aufgabe (c): 
    % ------------
    % Umwandlung in System 1. Ordnung:
    % --------------------------------
    %
    % y3 = y1' , y4 = y2'
    % 
    % => System 1. Ordnung:
    % y1' = y3                                        x1(0) = 0, 994
    % y2' = y4                                        x2(0) = 0
    % y3' = y1 + 2*y4 − mu'*(y1+mu)/D1 - mu*(y1-mu')/D2
    % y4' = y2 − 2*y3 − mu'*y2/D1-mu*y2/D2
    % 
    % Anfangswerte:
    % -------------
    % 
    % y1(0) = 0.994
    % y2(0) = 0
    % y3(0) = 0
    % y4(0) = −2.0015851063

    n=4;
    x0=0;
    xe=17.0652165602;
    mu = 0.012277471;
    muhat = 1 - mu; 


    m=10000;                                   
    h = (xe-x0)/m;                            
    y_euler = zeros(n,m+1);                       
    y_euler(:,1) = [0.994; 0; 0; -2.0015851063];   
    y_veuler = zeros(n,m+1);                       
    y_veuler(:,1) = [0.994; 0; 0; -2.0015851063]; 
    y_rk = zeros(n,m+1);                          
    y_rk(:,1) = [0.994; 0; 0; -2.0015851063];     

    x = x0;

    for it=1:m
        x=x0+it*h;                             
        y_euler(:,it+1) = euler(@fc,n,x,h,y_euler(:,it));
        y_veuler(:,it+1) = veuler(@fc,n,x,h,y_veuler(:,it));
        y_rk(:,it+1) = rungekutta(@fc,n,x,h,y_rk(:,it));
    end
    h1=h;
    lsg1_euler=y_euler;
    lsg1_veuler=y_veuler;
    lsg1_rk=y_rk;


    m=20000;                                  
    h = (xe-x0)/m; 
    y_euler = zeros(n,m+1);                       
    y_euler(:,1) = [0.994; 0; 0; -2.0015851063];   
    y_veuler = zeros(n,m+1);                       
    y_veuler(:,1) = [0.994; 0; 0; -2.0015851063]; 
    y_rk = zeros(n,m+1);                          
    y_rk(:,1) = [0.994; 0; 0; -2.0015851063];                                
    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(:,it+1) = euler(@fc,n,x,h,y_euler(:,it));
        y_veuler(:,it+1) = veuler(@fc,n,x,h,y_veuler(:,it));
        y_rk(:,it+1) = rungekutta(@fc,n,x,h,y_rk(:,it));
    end
    h2=h;
    lsg2_euler=y_euler;
    lsg2_veuler=y_veuler;
    lsg2_rk=y_rk;


    m=50000;                                
    h = (xe-x0)/m; 
    y_euler = zeros(n,m+1);                       
    y_euler(:,1) = [0.994; 0; 0; -2.0015851063];   
    y_veuler = zeros(n,m+1);                       
    y_veuler(:,1) = [0.994; 0; 0; -2.0015851063]; 
    y_rk = zeros(n,m+1);                          
    y_rk(:,1) = [0.994; 0; 0; -2.0015851063];     
    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(:,it+1) = euler(@fc,n,x,h,y_euler(:,it));
        y_veuler(:,it+1) = veuler(@fc,n,x,h,y_veuler(:,it));
        y_rk(:,it+1) = rungekutta(@fc,n,x,h,y_rk(:,it));
    end
    h3=h;
    lsg3_euler=y_euler;
    lsg3_veuler=y_veuler;
    lsg3_rk=y_rk;


    scrsz = get(0,'ScreenSize');
    figure('Position',[0 500 1400 500])
    subplot(1,3,1)
    plot(lsg1_euler(1,:),lsg1_euler(2,:),'k')
    hold on
    plot(lsg1_veuler(1,:),lsg1_veuler(2,:),'g')
    plot(lsg1_rk(1,:),lsg1_rk(2,:),'r')
    plot(1-mu,0,'*') % Mond
    plot(-mu,0,'o')  % Erde
    hold off
    xlabel('y_1(x)')
    ylabel('y_2(x)')
    axis square
    title('AWA 3, m= Schritte')
    subplot(1,3,2)
    plot(lsg2_euler(1,:),lsg2_euler(2,:),'k')
    hold on
    plot(lsg2_veuler(1,:),lsg2_veuler(2,:),'g')
    plot(lsg2_rk(1,:),lsg2_rk(2,:),'r')
    plot(1-mu,0,'*') % Mond
    plot(-mu,0,'o')  % Erde
    hold off
    xlabel('y_1(x)')
    ylabel('y_2(x)')
    axis square
    title('AWA 3, m= Schritte')
    subplot(1,3,3)
    plot(lsg3_euler(1,:),lsg3_euler(2,:),'k')
    hold on
    plot(lsg3_veuler(1,:),lsg3_veuler(2,:),'g')
    plot(lsg3_rk(1,:),lsg3_rk(2,:),'r')
    plot(1-mu,0,'*') % Mond
    plot(-mu,0,'o')  % Erde
    hold off
    xlabel('y_1(x)')
    ylabel('y_2(x)')
    axis square
    title('AWA 3, m= Schritte')
    legend('Euler','verbessertes Euler','Runge Kutta','Mond','Erde',2)


    % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %


    % Aufgabe (d): 
    % ------------

    n = 2; 
    x0 = 0; 
    xe = 30; 

    m=60;
    h = (xe-x0)/m; 
    y_euler = zeros(n,m+1);
    y_euler(1,1) = -1.03;
    y_euler(2,1) = 2.16;
    y_veuler = zeros(n,m+1);
    y_veuler(1,1) =  -1.03;
    y_veuler(2,1) = 2.16;
    y_rk = zeros(n,m+1);
    y_rk(1,1) =  -1.03;
    y_rk(2,1) = 2.16;

    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(:,it+1) = euler(@fd,n,x,h,y_euler(:,it));
        y_veuler(:,it+1) = veuler(@fd,n,x,h,y_veuler(:,it));
        y_rk(:,it+1) = rungekutta(@fd,n,x,h,y_rk(:,it));
    end

    h1=h;
    lsg1_euler=y_euler;
    lsg1_veuler=y_veuler;
    lsg1_rk=y_rk;

    scrsz = get(0,'ScreenSize');
    figure('Position',[0 500 1400 500])
    subplot(1,3,1)
    plot(x0:h1:xe,lsg1_euler(1,:),'k')
    hold on
    plot(x0:h1:xe,lsg1_euler(2,:),'k--')
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    legend('x(t)','y(t)',1)
    title('AWA 4, Euler Verfahren, x(t),y(t), m=60')
    subplot(1,3,2)
    plot(x0:h1:xe,lsg1_veuler(1,:),'g')
    hold on
    plot(x0:h1:xe,lsg1_veuler(2,:),'g--')
    hold off
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    legend('x(t)','y(t)',1)
    title('AWA 4, verbessteres Euler Verfahren, x(t),y(t), m=60')    
    subplot(1,3,3)
    plot(x0:h1:xe,lsg1_rk(1,:),'m')
    hold on
    plot(x0:h1:xe,lsg1_rk(2,:),'m--')
    hold off
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    legend('x(t)','y(t)',1)
    title('AWA 4, Runge Kutta Verfahren, x(t), y(t), m=60')
    
        
    n = 2; 
    x0 = 0; 
    xe = 30; 

    m=10000;
    h = (xe-x0)/m; 
    y_euler = zeros(n,m+1);
    y_euler(1,1) = -1.03;
    y_euler(2,1) = 2.16;
    y_veuler = zeros(n,m+1);
    y_veuler(1,1) =  -1.03;
    y_veuler(2,1) = 2.16;
    y_rk = zeros(n,m+1);
    y_rk(1,1) =  -1.03;
    y_rk(2,1) = 2.16;

    x = x0;

    for it=1:m
        x=x0+it*h;
        y_euler(:,it+1) = euler(@fd,n,x,h,y_euler(:,it));
        y_veuler(:,it+1) = veuler(@fd,n,x,h,y_veuler(:,it));
        y_rk(:,it+1) = rungekutta(@fd,n,x,h,y_rk(:,it));
    end

    h1=h;
    lsg1_euler=y_euler;
    lsg1_veuler=y_veuler;
    lsg1_rk=y_rk;

    scrsz = get(0,'ScreenSize');
    figure('Position',[0 500 1400 500])
    subplot(1,3,1)
    plot(x0:h1:xe,lsg1_euler(1,:),'k')
    hold on
    plot(x0:h1:xe,lsg1_euler(2,:),'k--')
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    legend('x(t)','y(t)',1)
    title('AWA 4, Euler Verfahren, x(t),y(t), m=10000')
    subplot(1,3,2)
    plot(x0:h1:xe,lsg1_veuler(1,:),'g')
    hold on
    plot(x0:h1:xe,lsg1_veuler(2,:),'g--')
    hold off
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    legend('x(t)','y(t)',1)
    title('AWA 4, verbessteres Euler Verfahren, x(t),y(t), m=10000')    
    subplot(1,3,3)
    plot(x0:h1:xe,lsg1_rk(1,:),'m')
    hold on
    plot(x0:h1:xe,lsg1_rk(2,:),'m--')
    hold off
    xlabel('t')
    ylabel('x(t),y(t)')
    axis square
    legend('x(t)','y(t)',1)
    title('AWA 4, Runge Kutta Verfahren, x(t), y(t), m=10000')

    
