function I = Gauss_quad(f, a, b, k)
%
% The function finds integral of f on the interval [a, b] using 
% Gaussian quadrature at k (k = 2, ..., 5) points  
%
t = [-0.5773502692  -0.7745966692   -0.8611363116   -0.9061798459;
    0.5773502692    0.0             -0.3399810436   -0.5384693101;
    0.0             0.7745966692    0.3399810436    0.0;
    0.0             0.0             0.8611363116    0.5384693101;
    0.0             0.0             0.0             0.9061798459]
c = [1.0            0.5555555556    0.3478548451    0.2369268850;
    1.0             0.8888888889    0.6521451549    0.4786286705;
    0.0             0.5555555556    0.6521451549    0.5688888889;
    0.0             0.0             0.3478548451    0.4786286705;
    0.0             0.0             0.0             0.2369268850]
%
% Transformation of the interval of integration
x(1 : k) = 0.5*((b-a).*t(1:k,k-1) + b + a);
%
y = feval(f, x);
%
cc(1 : k) = c(1 : k, k-1);
cd = cc';
%
int = y*cd;
I = int*(b-a)/2