function [Int] = Simpson_rule(fun,lim,n)
    % Integrates function fun using Simpson's rule
    % Input variables:
    % fun = integrand
    % lim = interval of integration
    % n     = number of nodes
    % Output variables:
    % Int = approximation to integral
    h = (lim(2)-lim(1))/n; x = linspace(lim(1),lim(2),n+1); y=fun(x);
    Int = h*(y(1)+4*sum(y(2:2:n))+2*sum(y(3:2:n-1))+y(n+1))/3;