Присмотритесь к формуле (1). Она говорит что интеграл можно представить как сумму членов, которые отличаются друг от друга только индексами переменной 
y. Причём индексы последовательно увеличиваются. 
Т.е. у нас должен быть цикл, в каждой итерации которого мы вычисляем 
yn и 
yn+1 и вычисляем член. Результат прибавляем к некоторой переменной. 
yn это ни что иное как результат вычисления функции 
f(x) в точке 
xn
Получается примерно так:
sum = 0;
h = (b-a)/1000;
for(x = a; x < b; x += h){
  sum += ( f(x) + f(x+h) )/2 * h;
}