function[y] = f(x)
    y = cos(x) / x;
end

function[wynik] = simpson(a, b, n)
    x = [a];
    y = [];
    i = 0;
    h = (b - a) / n;
    for i = 2:n+1
      x(i) = x(1) + (i - 1) * h;
    end

    for i = 1:n+1
      y(i) = f(x(i));
    end

    wynik = 0;

    for i = 1:n+1
      if i == 1 || i == n + 1
        wynik = wynik + y(i);
      else
        if mod(i, 2) == 0
          wynik = wynik + 4 * y(i);
        else
          wynik = wynik + 2 * y(i);
        end
      end
    end;

    wynik = wynik * (h / 3);
end