#include <iostream.h>
#include <math.h>

//==============================================================
void nrerror(char * error_text)
{

        cerr<<"Numerical Recipes run-time error : "<<error_text<<"...now exiting to system...\n";
        // exit(1);
}
//===========================================================================
double *vector(int N)

{
        double *v;
        v= new double [N];
        if (!v) nrerror("allocation failure in vector()");
 	return v;
}
//========================================
void free_vector(double *v)
{
        delete [] v;
}
//===========================================================================
void  rk4(double *y, double *dydx, int n, double x, double h, double *yout,
                void (*derivs)(double,double *,double *))
{
        int i;
        double xh,hh,h6,*dym,*dyt,*yt;
        dym=vector(n);
        dyt=vector(n);
        yt=vector(n);
        hh=h*0.5;
        h6=h/6.0;
        xh=x+hh;
        for (i=0;i<n;i++) yt[i]=y[i]+hh*dydx[i];
        (*derivs)(xh,yt,dyt);
        for (i=0;i<n;i++) yt[i]=y[i]+hh*dyt[i];
        (*derivs)(xh,yt,dym);
        for (i=0;i<n;i++) {
                yt[i]=y[i]+h*dym[i];
                dym[i] += dyt[i];
        }
        (*derivs)(x+h,yt,dyt);
        for (i=0;i<n;i++)
                yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
        free_vector(yt);
        free_vector(dyt);
        free_vector(dym);
}
//===========================================================================

#define PGROW -0.20
#define PSHRNK -0.25
#define FCOR 0.06666666         /* 1.0/15.0 */
#define SAFETY 0.9
#define ERRCON 6.0e-4

void  rkqc(double *y, double *dydx, int n, double *x, double htry, 
                double eps, double *yscal, double *hdid, double *hnext,
                void (*derivs)(double,double *,double *))
{
        int i;
        double xsav,hh,h,temp,errmax;
        double *dysav,*ysav,*ytemp;
        

        dysav=vector(n);
        ysav=vector(n);
        ytemp=vector(n);
        xsav=(*x);
        for (i=0;i<n;i++) {
                ysav[i]=y[i];
                dysav[i]=dydx[i];
        }
        h=htry;
        for (;;) {
                hh=0.5*h;
                rk4(ysav,dysav,n,xsav,hh,ytemp,derivs);
                *x=xsav+hh;
                (*derivs)(*x,ytemp,dydx);
                rk4(ytemp,dydx,n,*x,hh,y,derivs);
                *x=xsav+h;
                if (*x == xsav) nrerror("Step size too small in routine RKQC");
                rk4(ysav,dysav,n,xsav,h,ytemp,derivs);
                errmax=0.0;
                for (i=0;i<n;i++) {
                        ytemp[i]=y[i]-ytemp[i];
                        temp=fabs(ytemp[i]/yscal[i]);
                        if (errmax < temp) errmax=temp;
                }
                errmax /= eps;
                if (errmax <= 1.0) {
                        *hdid=h;
                        *hnext=(errmax > ERRCON ?
                                SAFETY*h*exp(PGROW*log(errmax)) : 4.0*h);
                        break;
                }
                h=SAFETY*h*exp(PSHRNK*log(errmax));
        }
        for (i=0;i<n;i++) y[i] += ytemp[i]*FCOR;
        free_vector(ytemp);
        free_vector(dysav);
        free_vector(ysav);
}

#undef PGROW
#undef PSHRNK
#undef FCOR
#undef SAFETY
#undef ERRCON
//===========================================================================

#define MAXSTP 10000
#define TINY 1.0e-30

int kmax=0,kount=0;  /* defining declaration */
double *xp=0,**yp=0,dxsav=0;  /* defining declaration */
int  odeint(double *ystart, int nvar, double x1, double x2, double eps,
                double h1, double hmin, int *nok, int *nbad,
                void (*derivs)(double,double *,double *))
{
        int nstp,i;
        double xsav,x,hnext,hdid,h;
        double *yscal,*y,*dydx;
       

        yscal=vector(nvar);
        y=vector(nvar);
        dydx=vector(nvar);
        x=x1;
        h=(x2 > x1) ? fabs(h1) : -fabs(h1);
        *nok = (*nbad) = kount = 0;
        for (i=0;i<nvar;i++) y[i]=ystart[i];
        if (kmax > 0) xsav=x-dxsav*2.0;
        for (nstp=1;nstp<=MAXSTP;nstp++) {
                (*derivs)(x,y,dydx);
                for (i=0;i<nvar;i++)
                        yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
                if (kmax > 0) {
                        if (fabs(x-xsav) > fabs(dxsav)) {
                                if (kount < kmax-1) {
                                        xp[++kount]=x;
                                        for (i=0;i<nvar;i++) yp[i][kount]=y[i];
                                        xsav=x;
                                }
                        }
                }
                if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;
                rkqc(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,derivs);
                if (hdid == h) ++(*nok); else ++(*nbad);
                if ((x-x2)*(x2-x1) >= 0.0) {
                        for (i=0;i<nvar;i++) ystart[i]=y[i];
                        if (kmax) {
                                xp[++kount]=x;
                                for (i=0;i<nvar;i++) yp[i][kount]=y[i];
                        }
                        free_vector(dydx);
                        free_vector(y);
                        free_vector(yscal);
                        return 0;
                }
                if (fabs(hnext) <= hmin) nrerror("Step size too small in ODEINT");
                h=hnext;
        }
        nrerror("Too many steps in routine ODEINT");
        return 1;
}
#undef MAXSTP
#undef TINY
//=================================================================

void RungeKutta(double *y, int N,double t1,double t2,void (*derivs)(double,double *,double *) ) 
{ 
        int nok,nbad; 
        odeint(y,N, t1, t2,1e-4,1e-2,0,&nok,&nbad,derivs); 
}