%
% Simpson Method for int(2xe^(-0.2x))dx;
% Finds the best h and n for given tolerance
%

f  = inline('2*x*exp(-0.2*x)');

n = 1;

%%%%%

a = 0;
b = 10;
actual = 50 - 150*exp(-2);

tol      = 10^(-5);
MAX_ITER = 40;

fprintf('------------------------------------------------\n');
fprintf('  n              h           sum          error \n');
fprintf('------------------------------------------------\n');
while( n < (MAX_ITER+1) )
    k = 2*n;
    h = (b-a)/k;
    sum = h*(f(a)+f(b))/3;
    if(n==1)
        x = a+h;
        sum = sum + (4*h/3)*f(x);
        dx = abs(sum-actual);
        fprintf('%3d   %8.6e   % 12.6f % 12.6e\n',n, h, sum, dx);
    else
        x = a+h;
        sum = sum + (4*h/3)*f(x);
        for j = 1:(n-1)
            x1 = a+(2*j)*h;
            x2 = a+(2*j+1)*h;
            sum = sum + (2*h/3)*f(x1)+(4*h/3)*f(x2);
        end
        dx = abs(sum-actual);
        fprintf('%3d   %8.6e   % 12.6f % 12.6e\n',n, h, sum, dx);
    end  
    if( dx<tol ) 
        return;
    end
    n = n+1;
end
fprintf('WARNING: Simpson method did not converge after %d steps.\n',n);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
