function integral = simpson( funktion, a, b, n )

    format long

    if mod( n, 2 ) ~= 0
        error( 'n needs to be an even number!' );
    end

    h = ( b - a ) / n;
    x = linspace( a, b, n + 1 );
    y = funktion( x );
    SumY=0;

    % Jämna värden
    for i = 1: 2: n - 2
           SumY = SumY + y( i ) + 4 * y( i + 1 ) + y( i + 2 );
    end

    integral = ( h / 3 ) * SumY;
end