function [x, y] = Trapezium(f, y0, a, b, h)
y(1) = y0;
n = (b - a)/h;
x = a:h:b;
for i=1:n
yt = y ( i ) + h * feval ( f, x ( i ), y ( i ) );
done = 0;
while ~done
y ( i + 1 ) = y ( i ) + 0.5 * h * ( feval ( f, x ( i ), y ( i ) )+feval ( f, x ( i+1 ), yt ) );
done = ( abs ( y ( i + 1 ) - yt ) < 1e-6 );
yt = y ( i + 1 );
end
end