function fun = myplots ( iexmpl, gdim, gscal, iprint )
% MyPlots
%
% This m-file is part of the material of the lectures on
%   Dynamical Systems and Nonlinear Analysis
% at the Math Department of Halle University, Germany.
%
% Author :      Martin Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Oct 27, 2011
%
% Wrapper for several plots  
%
% Reference:
%   A.M. Stuart, A.R. Humphries: Dynamical Systems and Numerical Analysis. -
%   Cambridge University Press - 1998.
%
% Parameters:
%   iexmpl (input)  : control flag example
%                       =  1 .. Example 2.1a), fixed point and oscillations with period 2
%                       =  2 .. Example 2.1b), quadratic map
%                       =  3 .. Example 2.1c), Henon map
%                       =  4 .. Remark 2.16(iii), invariants
%                       =  5 .. Remark 2.NN(i), stability of fixed points
%                       =  6 .. Example 2.NNb), stability of fixed points
%                       =  7 .. Example 2.NNb), stability of fixed points (ii)
%                       =  8 .. Example 2.NNc), stability of fixed points (iii)
%                       =  9 .. Remark 2.NN, strange attractors
%                       = 10 .. Remark 2.NN, strange attractors (ii)
%                       = 11 .. Example 2.14b), periodic solution
%                       = 12 .. Example 2.14c), aperiodic solution
%                       = 13 .. Example 2.14d), Lorenz attractor
%                       = 14 .. Example 2.14d), Lorenz attractor (ii)
%                       = 15 .. Example 2.14d), Lorenz attractor (iii)
%                       = 16 .. Remark 4.NN, time integration for hyperbolic stable periodic solution, heul = 0.10
%                       = 17 .. Remark 4.NN, time integration for hyperbolic stable periodic solution, heul = 0.20
%                       = 18 .. Remark 4.NN, time integration for hyperbolic stable periodic solution, heul = 0.35
%                       = 19 .. Remark 4.NN, time integration for hyperbolic stable periodic solution, heul = 0.36
%   gdim   (input)  : [ graphwidth graphheight] 
%                     dimension of graphical output in cm 
%   gscal  (input)  : [ facv shiftv fach shifth ] 
%                     rescaling factors for graphical output
%   iprint (input)  : control flag output
%                       = 0 .. screen output
%                       = 1 .. screen and file output
%
% Example:
%   myplots (  1, [ 24.0 10.0 ], [ 0.98  0.5 1.00 0.0 ], 0 );
%   myplots (  2, [ 24.0 10.0 ], [ 0.98  0.5 0.99 1.0 ], 0 );
%   myplots (  3, [ 24.0 10.0 ], [ 0.94  0.7 1.00 0.0 ], 0 );
%   myplots (  4, [ 24.0 10.0 ], [ 0.92  0.6 1.00 0.0 ], 0 );
%   myplots (  5, [ 24.0 10.0 ], [ 0.92  0.5 1.00 0.0 ], 0 );
%   myplots (  6, [ 24.0 10.0 ], [ 0.92  0.5 1.00 0.0 ], 0 );
%   myplots (  7, [ 24.0 10.0 ], [ 0.92  0.5 1.00 0.0 ], 0 );
%   myplots (  8, [ 24.0 10.0 ], [ 0.92  0.5 1.00 0.0 ], 0 );
%   myplots (  9, [ 24.0 10.0 ], [ 0.92  0.5 1.00 0.0 ], 0 );
%   myplots ( 10, [ 24.0 10.0 ], [ 0.92  0.5 1.00 0.0 ], 0 );
%   myplots ( 11, [ 24.0 10.0 ], [ 1.02 -0.5 1.00 0.0 ], 0 );
%   myplots ( 12, [ 24.0 10.0 ], [ 1.02 -1.0 1.00 0.0 ], 0 );
%   myplots ( 13, [ 24.0  9.0 ], [ 0.98  1.5 1.00 0.0 ], 0 );
%   myplots ( 14, [ 24.0  6.6 ], [ 0.92  1.4 1.00 0.0 ], 0 );
%   myplots ( 15, [ 24.0  6.6 ], [ 0.92  1.4 1.00 0.0 ], 0 );
%   myplots ( 16, [ 24.0 10.0 ], [ 0.95  0.5 1.00 0.0 ], 0 );

graphwidth  = gdim(1);
graphheight = gdim(2);

mygraph.facv    =  gscal(1);
mygraph.shiftv	=  gscal(2);
mygraph.fach	=  gscal(3);
mygraph.shifth	=  gscal(4);

papershift  =  0.0;

mygraph.ActFontSize      = 14;
mygraph.ActFontSizeSmall = 11;
mygraph.ActLineWidth     =  2;
mygraph.ActMarkerSize    =  6;
RatioFlag     =  0;

border      = 1.0;

%
% -> Bestimme Groesse des Bildschirms und Skalierungsfaktor fuer die Ausgabe
root = 0;
set ( root, 'Units', 'centimeters' );
screensize = get ( root, 'ScreenSize' );
left   = screensize(1);
bottom = screensize(2);
width  = screensize(3);
height = screensize(4);
set ( root, 'Units', 'pixels' );
if RatioFlag==0,
  screenscale = 1.0;
elseif RatioFlag==1,
  screenscale = (width-2*border) / graphwidth;
end;

%
% -> Definiere den Bildschirmausschnitt
if isempty ( get ( root, 'Children' ) ),
  h = figure;
else
  h = gcf;
  clf;
end;
set ( h, 'Units',    'centimeters' );
set ( h, 'Position', ...
  [ left+width-screenscale*graphwidth-border height-2*border-screenscale*graphheight ...
    screenscale*graphwidth                   screenscale*graphheight ] );
set ( h, 'Units',    'pixels' );

%
% -> Skalierung der Druckausgabe
set ( gcf, 'PaperUnits',    'centimeters' );
set ( gcf, 'PaperPosition', [ papershift papershift graphwidth graphheight ] );

% -> The plot
switch iexmpl
    case 1
        % -> Example 2.1a), fixed point and oscillations with period 2
        for iplot=1:2,
            switch iplot
                case 1
                    a    = 1.1;   U0   = 0.1;
                    Umin = 0.0;   Umax = 1.0;
                case 2
                    a    = 2.5;   U0   = 0.1;
                    Umin = 0.0;   Umax = 3.0;
            end
            N = 100;
            U = zeros ( N+1, 1);
            U(1) = U0;
            for n=1:N,
                U(n+1) = a * sin(U(n));
            end
            subplot ( 1, 2, iplot );
            plot ( 0:N, U, '.' );
            axis ( [ 0 N Umin Umax ] );
            xlabel ( 'n' );   ylabel ( 'U_n' );
            title ( sprintf ( 'a = %3.2f', a ) );
        end;
    case 2
        % -> Example 2.1b), quadratic map
        for iplot=1:2,
            switch iplot
                case 1
                    a    = 3.5;   U0   = 0.1;
                case 2
                    a    = 3.9;   U0   = 0.1;
            end
            Umin = 0.0;   Umax = 1.0;
            N = 100;
            U = zeros ( N+1, 1);
            U(1) = U0;
            for n=1:N,
                U(n+1) = a * U(n) * ( 1.0 - U(n) );
            end
            subplot ( 1, 2, iplot );
            plot ( 0:N, U, '.' );
            axis ( [ 0 N Umin Umax ] );
            xlabel ( 'n' );   ylabel ( 'U_n' );
            title ( sprintf ( 'Logistische Abbildung:  a = %0.3g', a ) );
        end
    case 3
        % -> Example 2.1c), Henon map
        a    =  0.0;   b    = 1.0;   X0   =  0.1;   Y0   = 0.1;
        Xmin = -0.4;   Xmax = 0.4;   Ymin = -0.4;   Ymax = 0.4;  
        tstr = sprintf ( 'Abbildung von Henon:  a = %0.3g, b = %0.3g', a, b );
        N = 500;
        X    = zeros ( N+1, 1);   Y    = zeros ( N+1, 1);
        X(1) = X0;                Y(1) = Y0;
        for n=1:N,
            X(n+1) = a - b*Y(n) - X(n)^2;
            Y(n+1) = X(n);
        end
        subplot ( 1, 2, 1 );
        plot ( X, Y, '.r' );
        axis ( [ Xmin Xmax Ymin Ymax ] );
        set ( gca, 'DataAspectRatio',    [ 1 1 1 ] );
        set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] );
        xlabel ( 'X_n' );   ylabel ( 'Y_n' );
        title ( tstr );
        subplot ( 1, 2, 2 );
        plot ( 0:N, X, '.', 'MarkerSize', 3 );
        axis ( [ 0 N Xmin Xmax ] );
        xlabel ( 'n' );   ylabel ( 'X_n' );
        title ( tstr );
    case 4
        % -> Remark 2.16(iii), invariants
        a    =  2.0;   X0   = 1.0;   Y0   =  1.0;
        Xmin = -2.5;   Xmax = 2.5;   Ymin = -2.5;   Ymax = 2.5;  
        tstr = sprintf('Invariante Menge  (x^2 + y^2)^{1/2} = %0.3g', a);
        N = 1000;
        X    = zeros ( N+1, 1);   Y    = zeros ( N+1, 1);
        X(1) = X0;                Y(1) = Y0;
        for n=1:N,
            r = sqrt(X(n)^2+Y(n)^2);
            X(n+1) = a * sin(r*X(n));
            Y(n+1) = r * cos(a*X(n));
        end
        subplot ( 1, 2, 1 );
        plot ( X, Y, '.r', 'MarkerSize', 4 );        
        hold on;
        for n=1:8,
            plot ( X(n), Y(n), '+r', 'MarkerSize', 9-n );
        end
        hold off;
        axis ( [ Xmin Xmax Ymin Ymax ] );
        set ( gca, 'DataAspectRatio',    [ 1 1 1 ] );
        set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] );
        xlabel ( 'X_n' );   ylabel ( 'Y_n' );
        title ( tstr );
        subplot ( 1, 2, 2 );
        plot ( 0:N, X, '.', 'MarkerSize', 2 );
        axis ( [ 0 N Xmin Xmax ] );
        xlabel ( 'n' );   ylabel ( 'X_n' );
        title ( tstr );
    case {5,6,7,8}
        % -> 5: Remark 2.NN(i), stability of fixed points
        % -> 6: Example 2.NNb), stability of fixed points
        % -> 7: Example 2.NNb), stability of fixed points (ii)
        % -> 8: Example 2.NNc), stability of fixed points (iii)
        for iplot=1:2,
            switch iexmpl
                case 5
                    switch iplot
                        case 1
                            X0   =  0.7;   Y0   = 0.8;
                            Xmin = -1.0;   Xmax = 1.0;   Ymin = -1.0;   Ymax = 1.0;
                        case 2
                            X0   =  0.75;  Y0   = 1.026;
                            Xmin = -10.0;   Xmax = 10.0;   Ymin = -10.0;   Ymax = 10.0;
                    end
                case {6,7}
                    switch iplot
                        case 1
                            X0   =  1.5;   Y0   = -1.5;
                            Xmin = -2.0;   Xmax =  2.0;   Ymin = -2.0;   Ymax = 2.0;
                        case 2
                            switch iexmpl
                                case 6
                                    X0   =  1.67;   Y0   = -1.67;
                                    Xmin = -10.0;   Xmax =  10.0;   Ymin = -10.0;   Ymax = 10.0;
                                case 7
                                    X0   =  0.82;   Y0   =  0.15;
                                    Xmin =  0.8;    Xmax =  1.2;    Ymin = -0.2;   Ymax = 0.2;
                            end
                    end
                case 8
                    switch iplot
                        case 1
                            X0   = -1.1;   Y0   = 0.1;
                            Xmin = -4.0;   Xmax = 0.0;   Ymin = -2.0;   Ymax = 2.0;
                        case 2
                            X0   = -1.001;   Y0   =  0.002;
                            Xmin = -1.5;     Xmax = -0.5;     Ymin = -0.5;   Ymax = 0.5;
                    end
            end
            N = 10;
            X    = zeros ( N+1, 1);   Y    = zeros ( N+1, 1);
            X(1) = X0;                Y(1) = Y0;
            switch iexmpl
                case 5
                    for n=1:N,
                        X(n+1) = 0.5*X(n) - Y(n)^3;
                        Y(n+1) = 0.5*Y(n) + X(n)*Y(n)^2;
                    end
                case {6,7}
                    for n=1:N,
                        X(n+1) = 1.25 - 0.25*X(n)^2 - Y(n)^2;
                        Y(n+1) = 0.5*X(n)*Y(n);
                    end
                case 8
                    for n=1:N,
                        X(n+1) = 1 - 2.0*X(n)^2 - Y(n)^2;
                        Y(n+1) = X(n)*Y(n);
                    end
            end
            subplot ( 1, 2, iplot );
            plot ( X(1), Y(1), '.r', 'MarkerSize', 16 );
            hold on;
            for n=2:N,
                plot ( X(n), Y(n), '+r', 'MarkerSize', N+1-n );
            end
            hold off;
            axis ( [ Xmin Xmax Ymin Ymax ] );
            set ( gca, 'DataAspectRatio',    [ 1 1 1 ] );
            set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] );
            xlabel ( 'X_n' );   ylabel ( 'Y_n' );
            title ( sprintf('X_0 = %0.4g,  Y_0 = %0.4g',X0,Y0) );
        end
    case {9,10}
        % ->  9: Remark 2.NN, strange attractors
        % -> 10: Remark 2.NN, strange attractors (ii)
        switch iexmpl
            case 9
                X0 = -0.5;   Y0 =  0.0;
            case 10
                X0 = -1.11;   Y0 = -0.408;
        end
        Xmin = -2.0;   Xmax =  2.0;   Ymin = -0.5;   Ymax = 0.5;  
        tstr = sprintf ( 'Seltsamer Attraktor:  X_0 = %0.3g, Y_0 = %0.3g', X0, Y0 );
        N = 5000;   N0 = 10;
        X    = zeros ( N+1, 1);   Y    = zeros ( N+1, 1);
        X(1) = X0;                Y(1) = Y0;
        for n=1:N,
            X(n+1) = 1 + Y(n) - 1.4*X(n)^2;
            Y(n+1) = 0.3*X(n);
        end
        subplot ( 1, 2, 1 );
        plot ( X(1), Y(1), '.r', 'MarkerSize', 16 );
        hold on;
        for n=2:N0,
            plot ( X(n), Y(n), '+r', 'MarkerSize', N0+1-n );
        end
        plot ( X((N0+1):(N+1)), Y((N0+1):(N+1)), '.m', 'MarkerSize', 1 );
        hold off,
        axis ( [ Xmin Xmax Ymin Ymax ] );
        set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] );
        xlabel ( 'X_n' );   ylabel ( 'Y_n' );
        title ( tstr );
        subplot ( 1, 2, 2 );
        plot ( 0:N0, X(1:(N0+1)), '+r', 'MarkerSize', 4 );
        hold on;
        plot ( (N0+1):(10*N0), X((N0+2):(10*N0+1)), '+m', 'MarkerSize', 4 );
        hold off;
        axis ( [ 0 10*N0 Xmin Xmax ] );
        xlabel ( 'n' );   ylabel ( 'X_n' );
        title ( tstr );
    case 11
        % -> 11: Example 2.14b), periodic solution
        X0    =  3.0;   Y0    = -3.0;
        tmin  = -0.0;   tmax  = 10.0;   
        xymin = -4.0;   xymax =  4.0;  
        nt = 1001;
        t = linspace ( tmin, tmax, nt );
        R = sqrt ( X0^2 + Y0^2 );
        r = exp(t)*R ./ sqrt ( 1 + R^2*(exp(2*t)-1) );
        x = r .* cos(-pi/4-t);
        y = r .* sin(-pi/4-t);
        subplot ( 1, 2, 1 );
        plot ( t, x, 'b', 'LineWidth', mygraph.ActLineWidth );
        hold on;
        plot ( t, y, 'g', 'LineWidth', mygraph.ActLineWidth );
        hold off,
        axis ( [ tmin tmax xymin xymax ] );
        xlabel ( 't' );   ylabel ( 'x(t), y(t)' );
        legend ( 'x(t)', 'y(t)' );
        subplot ( 1, 2, 2 );
        plot ( x, y, 'r', 'LineWidth', mygraph.ActLineWidth );
        axis ( [ xymin xymax xymin xymax ] );
        xlabel ( 'x(t)' );   ylabel ( 'y(t)' );
        set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] );
    case 12
        % -> 12: Example 2.14c), aperiodic solution
        tmin  = -0.0;   tmax  = 100.0;   
        xymin = -2.0;   xymax =   2.0;  
        umin  = -4.0;   umax  =   4.0;  
        nt = 1001;
        t = linspace ( tmin, tmax, nt );
        subplot ( 1, 2, 1 );
        plot ( exp(-t)+sin(t), exp(-t)+cos(10*pi/31*t), 'b'  ...
             , 'LineWidth', mygraph.ActLineWidth );
        axis ( [ xymin xymax xymin xymax ] );
        set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] );
        xlabel ( 'e^{-t}  +  sin t' );   
        ylabel ( 'e^{-t}  +  cos ( 10\pi * t / 31 )' );
        subplot ( 1, 2, 2 );
        plot ( t, exp(-t)+sin(t)+cos(sqrt(2)*t), 'r',   ...
             'LineWidth', mygraph.ActLineWidth );
        xlabel ( 't' );   ylabel ( 'e^{-t}  +  sin t  +  cos (2^{0.5} t)' );
        axis ( [ tmin tmax umin umax ] );
    case {13,14,15}
        % -> 13: Example 2.14d), Lorenz attractor
        % -> 14: Example 2.14d), Lorenz attractor (ii)
        % -> 15: Example 2.14d), Lorenz attractor (iii)
        tmin  = -0.0;   tmax  = 30.0;   
        X0 = 0.0;   Y0 = 10.0;   Z0 = 10.0;
        [ t, y ] = ode45 ( @mylorenz, [ tmin tmax ], [ X0, Y0, Z0 ]' );
        xmin = -20.0;   xmax =   20.0;  
        ymin = -30.0;   ymax =   30.0;  
        zmin =   0.0;   zmax =   50.0;  
        rmin = -10.0;   rmax =   40.0;  
        switch iexmpl
            case 13
                ThisLineWidth = mygraph.ActLineWidth;
                xxl    = t;      yyl    = y(:,1);   xxlstr = 't';    yylstr = 'x(t)';   
                xxlmin = tmin;   xxlmax = tmax;     yylmin = xmin;   yylmax = xmax;
                xxr    = t;      yyr    = y(:,3);   xxrstr = 't';    yyrstr = 'z(t)';
                xxrmin = tmin;   xxrmax = tmax;     yyrmin = zmin;   yyrmax = zmax;
            case 14
                ThisLineWidth = 1;
                xxl    = y(:,1);   yyl    = y(:,2);   xxlstr = 'x(t)';   yylstr = 'y(t)';   
                xxlmin = xmin;     xxlmax = xmax;     yylmin = ymin;     yylmax = ymax;
                xxr    = y(:,1);   yyr    = y(:,3);   xxrstr = 'x(t)';   yyrstr = 'z(t)';
                xxrmin = xmin;     xxrmax = xmax;     yyrmin = zmin;     yyrmax = zmax;
            case 15
                ThisLineWidth = 1;
                xxl    = y(:,2);   yyl    = y(:,3);   xxlstr = 'y(t)';   yylstr = 'z(t)';   
                xxlmin = ymin;     xxlmax = ymax;     yylmin = zmin;     yylmax = zmax;
                xxr    = y(:,3);   yyr    = sqrt(y(:,1).^2+y(:,2).^2);   
                xxrstr = 'z(t)';   yyrstr = '( x^2(t) + y^2(t) )^{0.5}';
                xxrmin = zmin;     xxrmax = zmax;     yyrmin = rmin;     yyrmax = rmax;
        end
        subplot ( 1, 2, 1 );
        plot ( xxl, yyl, 'LineWidth', ThisLineWidth );
        axis ( [ xxlmin xxlmax yylmin yylmax ] );
        xlabel ( xxlstr );   
        ylabel ( yylstr );   
        subplot ( 1, 2, 2 );
        plot ( xxr, yyr, 'LineWidth', ThisLineWidth );
        axis ( [ xxrmin xxrmax yyrmin yyrmax ] );
        xlabel ( xxrstr );   
        ylabel ( yyrstr );   
    case {16,17,18,19}
        % -> {16,17,18,19} .. Remark 4.NN, time integration for hyperbolic stable periodic solution
        %     16           .. heul = 0.10
        %     17           .. heul = 0.20
        %     18           .. heul = 0.35
        %     19           .. heul = 0.36
        switch iexmpl
            case 16
                heul = 1.0e-1;
            case 17
                heul = 2.0e-1;
            case 18
                heul = 3.5e-1;
            case 19
                heul = 3.6e-1;
        end
        X0    =  1.8;   Y0    = -1.8;
        tmin  = -0.0;   tmax  = 20.0;   
        xymin = -2.0;   xymax =  2.0;  
        nfirst = 1:5;
        MarkerSizeFirst = 8;
        atoleulim = 1.0e-12;
        rtoleulim = 1.0e-10;
        % -> analytical solution
        nt = 1001;
        t = linspace ( tmin, tmax, nt );
        R = sqrt ( X0^2 + Y0^2 );
        r = exp(t)*R ./ sqrt ( 1 + R^2*(exp(2*t)-1) );
        x = r .* cos(-pi/4-t);
        y = r .* sin(-pi/4-t);
        % -> numerical solution
        teul  = tmin:heul:tmax;
        Neul  = length(teul) - 1;
        % -> Explicit Euler method
        ycur  = [ X0 Y0 ]';
        yeulex = zeros ( 1+Neul, 2 );
        yeulex(1,:) = ycur';
        for ieul=1:Neul
            tcur = teul(ieul);
            fcur = myrhsper ( tcur, ycur );
            ycur = ycur + heul * fcur;
            yeulex(1+ieul,:) = ycur';
        end
%         % -> Implicit Euler method
%         ycur  = [ X0 Y0 ]';
%         yeulim = zeros ( 1+Neul, 2 );
%         yeulim(1,:) = ycur';
%         for ieul=1:Neul
%             tcur = teul(ieul);
%             fcur = myrhsper ( tcur, ycur );
%             yit  = ycur + heul * fcur;
%             diff = norm(heul*fcur);
%             while diff>=atoleulim+rtoleulim*abs(ycur),
%                 fcur = myrhsper ( tcur, yit );
%                 diff = norm ( ycur + heul * fcur - yit );
%                 yit  = ycur + heul * fcur;
%             end
%             ycur = yit;
%             yeulim(1+ieul,:) = ycur';
%         end
        % -> phase plot
        hphaseplot = subplot ( 1, 2, 1 );
        hphaseexact = plot ( x, y, 'b', 'LineWidth', mygraph.ActLineWidth+1 );
        hold on
        hphaseeulex = plot ( yeulex(:,1), yeulex(:,2), 'r', 'LineWidth', mygraph.ActLineWidth );        
        plot ( yeulex(nfirst,1), yeulex(nfirst,2), 'xr', 'MarkerSize', MarkerSizeFirst );
%         hphaseeulim = plot ( yeulim(:,1), yeulim(:,2), 'g', 'LineWidth', mygraph.ActLineWidth );        
%         plot ( yeulim(nfirst,1), yeulim(nfirst,2), 'xg', 'MarkerSize', MarkerSizeFirst );
        hold off
        axis ( [ xymin xymax xymin xymax ] );
        xlabel ( 'x(t)' );   ylabel ( 'y(t)' );
        title ( sprintf('Explicit Euler method,  h = %3.2f',heul) );
        set ( gca, 'PlotBoxAspectRatio', [ 1 1 1 ] );
        % -> plot of time history
        htimeplot = subplot ( 1, 2, 2 );
        htimeexact = plot ( t, x, 'b', 'LineWidth', mygraph.ActLineWidth );
        hold on;
        htimeeulex = plot ( teul, yeulex(:,1), 'r', 'LineWidth', mygraph.ActLineWidth );        
        plot ( teul(nfirst), yeulex(nfirst,1), 'xr', 'MarkerSize', MarkerSizeFirst );
%         htimeeulim = plot ( teul, yeulim(:,1), 'g', 'LineWidth', mygraph.ActLineWidth );        
%         plot ( teul(nfirst), yeulim(nfirst,1), 'xg', 'MarkerSize', MarkerSizeFirst );
        hold off,
        axis ( [ tmin tmax xymin xymax ] );
        xlabel ( 't' );   ylabel ( 'x(t)' );
        title ( sprintf('Explicit Euler method,  h = %3.2f',heul) );
end

%
% -> Korrektur der Groesse und Lage des Koordinatensystems
allch = get ( gcf, 'Children' );
for ax=allch',
    rect    = get ( ax, 'Position' );
    rect(2) = rect(2) + mygraph.shiftv * (1-mygraph.facv) * rect(4);
    rect(4) = mygraph.facv * rect(4);
    rect(1) = rect(1) + mygraph.shifth * (1-mygraph.fach) * rect(3);
    rect(3) = mygraph.fach * rect(3);
    set ( ax, 'Position', rect );
    set ( ax, 'FontSize', mygraph.ActFontSize );
    set ( get ( ax, 'XLabel' ), 'FontSize', mygraph.ActFontSize );
    set ( get ( ax, 'YLabel' ), 'FontSize', mygraph.ActFontSize );
    set ( get ( ax, 'Title' ), 'FontSize', mygraph.ActFontSize );
end

switch iexmpl
    case {16,17,18,19}
        for allax=[hphaseplot htimeplot],
            set ( gcf, 'CurrentAxes', allax );
            set ( gca, 'FontSize', mygraph.ActFontSizeSmall );
            if (allax==hphaseplot),
                hleg = legend ( [ hphaseexact hphaseeulex ] ...
                              , 'Analytical solution'       ...
                              , 'Numerical solution'        ...
                              , 'Location', 'NorthEast'     ...
                              );
            elseif (allax==htimeplot),
                hleg = legend ( [ htimeexact htimeeulex ] ...
                              , 'Analytical solution'     ...
                              , 'Numerical solution'      ...
                              , 'Location', 'SouthEast'   ...
                              );
            end

            rect = get ( hleg, 'Position' );
            set ( gca, 'FontSize', mygraph.ActFontSize );

            allch = get ( hleg, 'Children' );
            for al=allch',
                if strcmp(get(al,'Type'),'text'),
                    set ( al, 'FontSize', mygraph.ActFontSizeSmall );
                end;
            end;
            % rect(1) = rect(1) - mygraph.shifthleg*rect(3);
            set ( hleg, 'Position', rect );
            set ( hleg, 'Position', rect );

        end
end

%
% -> output to file (if applicable)
if iprint==1
    if iexmpl<10,
        iexstr = sprintf ( '0%1i', iexmpl );
    else
        iexstr = sprintf ( '%2i',  iexmpl );
    end
    eval ( [ 'print -dpng ../png/myplots' iexstr '.png' ] );
end

%------------------------------------------------------------------------------

function yprime = mylorenz ( t, y )

sigma = 10.0;
r     = 28.0;
b     = 8.0/3.0;

yprime = zeros ( size ( y ) );

yprime(1) = sigma * ( y(2) - y(1) );
yprime(2) = r*y(1) - y(2) - y(1)*y(3);
yprime(3) = y(1)*y(2) - b*y(3);

%------------------------------------------------------------------------------

function yprime = myrhsper ( t, y )

yprime = zeros ( size ( y ) );

u = y(1);
v = y(2);

yprime(1) =  u + v - u*(u^2+v^2);
yprime(2) = -u + v - v*(u^2+v^2);

