// ---------------------------------------------------------
// Fonction integration trapezes
// ---------------------------------------------------------
function [v]=integtrap(f,a,b,N); 
//
h = (b-a)/N;
pts = linspace(a,b,N+1); // = a + h * (0:N);
v = h * ( sum(f(pts(2:$-1))) + 0.5 * (f(pts(1)) + f(pts($))) );
//
endfunction
//
// ---------------------------------------------------------
// Fonction integration trapezes ou Simpson
// ---------------------------------------------------------
function [v]=integsymp(f,a,b,N); 
//
if (modulo(N,2)~=0) then
 N=N+1;
end
h = (b-a)/N;
pts = linspace(a,b,N+1); // = a + h * (0:N);
v = (h/3)*( f(pts(1)) + f(pts($)) + 4*sum(f(pts(2:2:$-1))) + 2*sum(f(pts(3:2:$-2))));
//
endfunction
//
// ---------------------------------------------------------
// Fonctions tests 
// ---------------------------------------------------------
//
function [y]=f1x2(x);
y = (1+x.*x).^(-1);
endfunction
//
function [y]=f1sin2x(x);
y = (1+sin(x).*sin(x)).^(-1);
endfunction
//
function [y]=fex1x2(x);
y=exp(-x).*((1+x.*x).^(-1));
endfunction
//