void qsimp(double a, double b, double& s)
{
    double eps = 0.000001;
    const int jmax = 20;
    double ost = -1e+30;
    double st, os = -1e+30;
    for (int j = 1; j<=jmax; j++)
	{
        trapzd(a, b, st, j);
        s = (4.0 * st - ost) / 3.0;
        if (fabs(s - os) < eps * fabs(os))
		{
			return;
		}
        os = s;
        ost = st;
    }
    cout<<"too many steps."<<endl;
}