function I = trap_int(func,a,b,N)
    
    h = (b-a)/N;
    x = a:h:b;
    
    f = func(x);
    
    I = (f(1)+f(end))*h/2 + h*sum(f(2:end-1));
    