Метод трапеций

#include <iostream>
#include <cstdlib>
 
using namespace std;
const double PRECISION = 0.001; // точность вычисления
 
/* фунция вычисления площади трапеции */
inline double getTrapezoidArea( double a, double b, double h );
 
int main()
{
   double a = 0.0, b = 0.0, c = 0.0, d = 0.0, e = 0.0; // коэффициенты
   double y1 = 0.0, y2 = 0.0;
   double downLimit = 0.0, upLimit = 0.0; // пределы интегрирования
   double S = 0.0; // площадь под графиком
 
   cout << "Добро пожаловать в программу вычисления интеграла полинома 4-ой степени от a до b методом трапеции с точностью " << PRECISION << endl;
   cout << "Введите коэффициенты полинома ax^4 + bx^3 +cx^2 +dx + e" << endl;
 
   cout << "a = ";
   cin >> a;
   if ( a == 0.0 )
   {
      cout << "Ошибка. Это не полином четвертой степени. Выход" << endl;
      exit( EXIT_SUCCESS );
   }
   cout << "b = ";
   cin >> b;
   cout << "c = ";
   cin >> c;
   cout << "d = ";
   cin >> d;
   cout << "e = ";
   cin >> e;
 
   cout << "Введите пределы интегрирования( от a до b, где a <= b )" << endl;
   cout << "Введите нижний предел: ";
   cin >> downLimit;
   cout << "Введите верхний предел: ";
   cin >> upLimit;
 
   if ( downLimit > upLimit )
   {
      cout << "Ошибка. Неправильно заданы пределы интегрирования. Выход" << endl;
      exit( EXIT_SUCCESS );
   }
   else if ( downLimit == upLimit )
   {
      S = 0.0;
   }
   else
   {
      S = 0.0;
      y1 = a * downLimit*downLimit*downLimit*downLimit
           + b * downLimit*downLimit*downLimit
           + c * downLimit*downLimit + d * downLimit + e;
 
      for ( double i = downLimit + PRECISION; i <= upLimit; i += PRECISION )
      {
         y2 = a * i*i*i*i + b * i*i*i + c * i*i + d * i + e;
         if ( y1 >= 0.0 && y2 >= 0.0 )
         {
            S += getTrapezoidArea( y1, y2, PRECISION );
         }
         y1 = y2;
      }
   }
 
   cout << "Результат: " << S << endl;
   return 0;
}
 
inline double getTrapezoidArea( double a, double b, double h )
{
   return ( a + b ) * h / 2.0;
}

S=0.5*(func(upLimit)+func(downLimit));
for ( double i = downLimit + PRECISION; i <= upLimit - PRECISION ; i += PRECISION )
            S += func(i);
 S=S*PRECISION;