function plot_waves
% long transient metastable oscillations
% 2013. 12. 10.

    format long;
    
    
    f = @(x) 0.5.*(abs(x+1) - abs(x-1));
    
    
    alfa = 3.5;
    beta = 2.5;
    N_arraysize = 8; % 6-8-16

    %8:
    % A =[b c 0 0 0 0 0 a; ...
    %     a b c 0 0 0 0 0; ...
    %     0 a b c 0 0 0 0; ...
    %     0 0 a b c 0 0 0; ...
    %     0 0 0 a b c 0 0; ...
    %     0 0 0 0 a b c 0; ...
    %     0 0 0 0 0 a b c; ...
    %     c 0 0 0 0 0 a b];

    a = alfa;
    b = 0.0;
    c = beta;
    
    multiplier_matrix = zeros(N_arraysize);
    for p=1:N_arraysize
        for q=1:N_arraysize
            if (p==q)
                multiplier_matrix(p, q) = b;
            end
            if (p==q-1)
                multiplier_matrix(p, q) = c;
            end
            if (p==q+1)
                multiplier_matrix(p, q) = a;
            end
            if (p==1 && q==N_arraysize)
                multiplier_matrix(p, q) = a;
            end
            if (p==N_arraysize && q==1)
                multiplier_matrix(p, q) = c;
            end
        end
    end
    
    
    simulation_length = 50;

    y0 = [];
    
    %{
    init_state = '{N/2}+{N/2}-';
    tmp = 1.0;
    for i=1:floor(N_arraysize/2)
        y0 = [y0; tmp];
    end
    for i=1:ceil(N_arraysize/2)
        y0 = [y0; -tmp];
    end
    %}
    
    init_state = '5+3-';
    y0 = [1 1 1 1 1 -1 -1 -1];


    %%{
    %%% explicit Euler
    
    %%}
    
    %
    %%% ode45
    options = odeset('RelTol', 1e-7, 'AbsTol', 1e-20);
    tt = [0, simulation_length];
    [T Y] = ode45(@(t, y) state_equations(y, f, multiplier_matrix), tt, y0, options);
    %}


    figure;
    hold on;
    grid on;
    plot(T, Y(:, 1), 'r.');
    plot(T, Y(:, 2), 'g.');
    plot(T, Y(:, 3), 'b.');
    plot(T, Y(:, 4), 'm.');
    plot(T, Y(:, 5), 'y.');
    plot(T, Y(:, 6), 'c.');
    plot(T, Y(:, 7), 'k.');
    plot(T, Y(:, 8), '.', 'Color', [0.7 0 0]);
    %plot(T, Y(:, 9), '.', 'Color', [1 0.8 0]);
    %plot(T, Y(:, 10), '.', 'Color', [0.7 0.3 0.7]);
    plot([T(1) T(end)], [1 1], 'k-');
    plot([T(1) T(end)], [-1 -1], 'k-');
    title(strcat('alfa=', num2str(alfa), ', beta=', num2str(beta), ', cell-num=', num2str(N_arraysize), ', init-state=', init_state));
    set(gca, 'YLim', [-7 7]);
    hold off;





end
