function main
    disp '---------------------------------';

    function [ok] = myplot(num1, den1, desc1, num2, den2, desc2)
        figure;

        subplot(1, 2, 1);
        step(num1, den1);
        hold on;
        step(num2, den2);
        hold off;
        legend(desc1, desc2);

        subplot(1, 2, 2);
        impulse(num1, den1);
        hold on;
        impulse(num2, den2);
        hold off;
        legend(desc1, desc2);
    end

    % inercyjny I-go rzedu
    k1 = 1; T1 = 2;
    k2 = 3; T2 = 1;
    num1 = [0, k1]; den1 = [T1, 1]; desc1 = ['k = ', num2str(k1), ', T = ', num2str(T1)];
    num2 = [0, k2]; den2 = [T2, 1]; desc2 = ['k = ', num2str(k2), ', T = ', num2str(T2)];
    myplot(num1, den1, desc1, num2, den2, desc2);

    % inercyjny II-go rzedu
    k1 = 1; T11 = 2; T21 = 3;
    k2 = 2; T12 = 1; T22 = 5;
    num1 = [0, 0, k1]; den1 = [T11 * T21, T11 + T21, 1]; desc1 = ['k = ', num2str(k1), ', T1 = ', num2str(T11), ', T2 = ', num2str(T21)];
    num2 = [0, 0, k2]; den2 = [T12 * T22, T12 + T22, 1]; desc2 = ['k = ', num2str(k2), ', T1 = ', num2str(T12), ', T2 = ', num2str(T22)];
    myplot(num1, den1, desc1, num2, den2, desc2);

    % inercyjny II-go rzedu (oscylacyjny)
    k1 = 1; T1 = 2; xi1 = 0.5;
    k2 = 2; T2 = 1; xi2 = 1.5;
    num1 = [0, 0, k1]; den1 = [T1^2, 2*xi1*T1, 1]; desc1 = ['k = ', num2str(k1), ', T = ', num2str(T1), ', \xi = ', num2str(xi1)];
    num2 = [0, 0, k2]; den2 = [T2^2, 2*xi2*T2, 1]; desc2 = ['k = ', num2str(k2), ', T = ', num2str(T2), ', \xi = ', num2str(xi2)];
    myplot(num1, den1, desc1, num2, den2, desc2);

    % calkujacy idealny
    k1 = 1; Ti1 = 2;
    k2 = 2; Ti2 = 1;
    num1 = [0, k1]; den1 = [Ti1, 0]; desc1 = ['k = ', num2str(k1), ', Ti = ', num2str(Ti1)];
    num2 = [0, k2]; den2 = [Ti2, 0]; desc2 = ['k = ', num2str(k2), ', Ti = ', num2str(Ti2)];
    myplot(num1, den1, desc1, num2, den2, desc2);

    % calkujacy rzeczywisty
    k1 = 1; Ti1 = 4; T1 = 2;
    k2 = 2; Ti2 = 5; T2 = 3;
    num1 = [0, 0, k1]; den1 = [Ti1 * T1, Ti1, 0]; desc1 = ['k = ', num2str(k1), ', Ti = ', num2str(Ti1), ', T = ', num2str(T1)];
    num2 = [0, 0, k2]; den2 = [Ti2 * T2, Ti2, 0]; desc2 = ['k = ', num2str(k2), ', Ti = ', num2str(Ti2), ', T = ', num2str(T2)];
    myplot(num1, den1, desc1, num2, den2, desc2);

    % różniczkujący idealny
    %
    % niemożliwy w ten sposób

    % różniczkujący rzeczywisty
    Td1 = 2; T1 = 0.5;
    Td2 = 3; T2 = 0.25;
    num1 = [Td1, 0]; den1 = [T1, 1]; desc1 = ['Td = ', num2str(Td1), ', T = ', num2str(T1)];
    num2 = [Td2, 0]; den2 = [T2, 1]; desc2 = ['Td = ', num2str(Td2), ', T = ', num2str(T2)];
    myplot(num1, den1, desc1, num2, den2, desc2);

    % inercyjny I-go rzedu z opoznieniem
    k1 = 3; T1 = 0.5; theta1 = 2; n1 = 5;
    k2 = 2; T2 = 2;   theta2 = 5; n2 = 7;
    num_i1 = [0, k1]; den_i1 = [T1, 1]; desc1 = ['k = ', num2str(k1), ', T = ', num2str(T1), ', \theta = ', num2str(theta1), ', n = ', num2str(n1)];
    num_i2 = [0, k2]; den_i2 = [T2, 1]; desc2 = ['k = ', num2str(k2), ', T = ', num2str(T2), ', \theta = ', num2str(theta2), ', n = ', num2str(n2)];
    [num_d1, den_d1] = pade(theta1, n1);
    [num_d2, den_d2] = pade(theta2, n2);
    [num1, den1] = series(num_d1, den_d1, num_i1, den_i1);
    [num2, den2] = series(num_d2, den_d2, num_i2, den_i2);
    myplot(num1, den1, desc1, num2, den2, desc2);
end