int m,i,j,k;
    double hh,p,dt,x,tt,q,a[4],*g,*b,*c,*d,*e;
    g=(double*)malloc(n*sizeof(double));
    b=(double*)malloc(n*sizeof(double));
    c=(double*)malloc(n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    e=(double*)malloc(n*sizeof(double));
    hh=h; m=1; p=1.0+eps; x=t;
    for (i=0; i<=n-1; i++) c[i]=y[i];
    while (p>=eps)
      { a[0]=hh/2.0; a[1]=a[0]; a[2]=hh; a[3]=hh;
        for (i=0; i<=n-1; i++)
          { g[i]=y[i]; y[i]=c[i];}
        dt=h/m; t=x;
        for (j=0; j<=m-1; j++)
          { ghamgf(t,y,n,d);
            for (i=0; i<=n-1; i++) 
              { b[i]=y[i]; e[i]=y[i];}
            for (k=0; k<=2; k++)
              { for (i=0; i<=n-1; i++)
                  { y[i]=e[i]+a[k]*d[i];
                    b[i]=b[i]+a[k+1]*d[i]/3.0;
                  }
                tt=t+a[k];
                ghamgf(tt,y,n,d);
              }
            for (i=0; i<=n-1; i++)
              y[i]=b[i]+hh*d[i]/6.0;
            t=t+dt;
          }
        p=0.0;
        for (i=0; i<=n-1; i++)
          { q=fabs(y[i]-g[i]);
            if (q>p) p=q;
          }
        hh=hh/2.0; m=m+m;
      }
    free(g); free(b); free(c); free(d); free(e);
    return;
}
void ghamg(double t,double h,int n,double y[],double eps,int k,double z[],
		void (*ghamgf)(double t,double y[],int n,double d[]))
{ 
	//extern void ghamgf();
    //void rkt5();
    int i,j,m;
    double a,q,*b,*d,*u,*v,*w,*g;
    b=(double*)malloc(4*n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    u=(double*)malloc(n*sizeof(double));
    v=(double*)malloc(n*sizeof(double));
    w=(double*)malloc(n*sizeof(double));
    g=(double*)malloc(n*sizeof(double));
    a=t;
    for (i=0; i<=n-1; i++) z[i*k]=y[i];
    ghamgf(t,y,n,d);
    for (i=0; i<=n-1; i++) b[i]=d[i];
    for (i=1; i<=3; i++)
      if (i<=k-1)
        { t=a+i*h;
          rkt5(t,h,y,n,eps,ghamgf);
          for (m=0; m<=n-1; m++) z[m*k+i]=y[m];
          ghamgf(t,y,n,d);
          for (m=0; m<=n-1; m++) b[i*n+m]=d[m];
        }
    for (i=0; i<=n-1; i++) u[i]=0.0;
    for (i=4; i<=k-1; i++)
      { for (j=0; j<=n-1; j++)
          { q=2.0*b[3*n+j]-b[n+n+j]+2.0*b[n+j];
            y[j]=z[j*k+i-4]+4.0*h*q/3.0;
          }
        for (j=0; j<=n-1; j++)
          y[j]=y[j]+112.0*u[j]/121.0;
        t=a+i*h;
        ghamgf(t,y,n,d);
        for (j=0; j<=n-1; j++)
          { q=9.0*z[j*k+i-1]-z[j*k+i-3];
            q=(q+3.0*h*(d[j]+2.0*b[3*n+j]-b[n+n+j]))/8.0;
            u[j]=q-y[j];
            z[j*k+i]=q-9.0*u[j]/121.0;
            y[j]=z[j*k+i];
            b[n+j]=b[n+n+j];
            b[n+n+j]=b[n+n+n+j];
          }
        ghamgf(t,y,n,d);
        for (m=0; m<=n-1; m++) b[3*n+m]=d[m];
      }
    free(b); free(d); free(u); free(v); free(w); free(g);
    return;
}
/////////////////////////////////////////////////////////////
void gtnr1(double t,double h,int n,double y[],
		void (*gtnr1f)(double t,double y[],int n,double d[]))
{ 
	//extern void gtnr1f();
    int j;
    double s,aa,bb,dd,g,dy,dy1,*d,*p,*w,*q,*r;
    w=(double*)malloc(4*n*sizeof(double));
    q=(double*)malloc(4*n*sizeof(double));
    r=(double*)malloc(4*n*sizeof(double));
    d=(double*)malloc(n*sizeof(double));
    p=(double*)malloc(n*sizeof(double));
    for (j=0; j<=n-1; j++) w[j]=y[j];
    gtnr1f(t,y,n,d);
    for (j=0; j<=n-1; j++)
      { q[j]=d[j]; y[j]=w[j]+h*d[j]/2.0;
        w[n+j]=y[j];
      }
    s=t+h/2.0;
    gtnr1f(s,y,n,d);
    for (j=0; j<=n-1; j++)
      { q[n+j]=d[j];
        y[j]=w[j]+h*d[j]/2.0;
        w[n+n+j]=y[j];
      }
    gtnr1f(s,y,n,d);
    for (j=0; j<=n-1; j++) q[n+n+j]=d[j];
    for (j=0; j<=n-1; j++)
      { aa=q[n+n+j]-q[n+j];
        bb=w[n+n+j]-w[n+j];
        if (-aa*bb*h>0.0)
          { p[j]=-aa/bb; dd=-p[j]*h;
            r[j]=exp(dd);
            r[n+j]=(r[j]-1.0)/dd;
            r[n+n+j]=(r[n+j]-1.0)/dd;
            r[3*n+j]=(r[n+n+j]-1.0)/dd;
          }
        else p[j]=0.0;
        if (p[j]<=0.0) g=q[n+n+j];
        else
          { g=2.0*(q[n+n+j]-q[j])*r[n+n+j];
            g=g+(q[j]-q[n+j])*r[n+j]+q[n+j];
          }
        w[3*n+j]=w[j]+g*h;
        y[j]=w[3*n+j];
      }
    s=t+h;
    gtnr1f(s,y,n,d);
    for (j=0; j<=n-1; j++) q[3*n+j]=d[j];
    for (j=0; j<=n-1; j++)
      { if (p[j]<=0.0)
          { dy=q[j]+2.0*(q[n+j]+q[n+n+j]);
            dy=(dy+q[n+n+n+j])*h/6.0;
          }
        else
          { dy=-3.0*(q[j]+p[j]*w[j])+2.0*(q[n+j]
                          +p[j]*w[n+j]);
            dy=dy+2.0*(q[n+n+j]+p[j]*w[n+n+j]);
            dy=dy-(q[n+n+n+j]+p[j]*w[n+n+n+j]);
            dy=dy*r[n+n+j]+q[j]*r[n+j];
            dy1=q[j]-q[n+j]-q[n+n+j]+q[n+n+n+j];
            dy1=dy1+(w[j]-w[n+j]-w[n+n+j]+w[n+n+n+j])*p[j];
            dy=(dy+4.0*dy1*r[n+n+n+j])*h;
          }
        y[j]=w[j]+dy;
      }
    free(d); free(p); free(w); free(q); free(r);
    return;
}
/////////////////////////////////////////////////////////////
void gtnr2(double t,double h,int n,double y[],int k,double z[],
		void (*gtnr2f)(double t,double y[],int n,double d[]))
{ 
	//extern void gtnr2f();
    int i,j;
    double a,s,aa,bb,dd,g,dy,dy1,*d,*p,*w,*q,*r;
    d=(double*)malloc(n*sizeof(double));
    p=(double*)malloc(n*sizeof(double));
    w=(double*)malloc(4*n*sizeof(double));
    q=(double*)malloc(4*n*sizeof(double));
    r=(double*)malloc(4*n*sizeof(double));
    a=t;
    for (j=0; j<=n-1; j++) z[j*k]=y[j];
    for (i=1; i<=k-1; i++)
      { t=a+(i-1)*h;
        for (j=0; j<=n-1; j++) w[j]=y[j];
        gtnr2f(t,y,n,d);
        for (j=0; j<=n-1; j++)
          { q[j]=d[j]; y[j]=w[j]+h*d[j]/2.0;
            w[n+j]=y[j];
          }
        s=t+h/2.0;
        gtnr2f(s,y,n,d);
        for (j=0; j<=n-1; j++)
          { q[n+j]=d[j];
            y[j]=w[j]+h*d[j]/2.0;
            w[n+n+j]=y[j];
          }
        gtnr2f(s,y,n,d);
        for (j=0; j<=n-1; j++) q[n+n+j]=d[j];
        for (j=0; j<=n-1; j++)
          { aa=q[n+n+j]-q[n+j];
            bb=w[n+n+j]-w[n+j];
            if (-aa*bb*h>0.0)
              { p[j]=-aa/bb; dd=-p[j]*h;
                r[j]=exp(dd);
                r[n+j]=(r[j]-1.0)/dd;
                r[n+n+j]=(r[n+j]-1.0)/dd;
                r[3*n+j]=(r[n+n+j]-1.0)/dd;
              }
            else p[j]=0.0;
            if (p[j]<=0.0) g=q[n+n+j];
            else
              { g=2.0*(q[n+n+j]-q[j])*r[n+n+j];
                g=g+(q[j]-q[n+j])*r[n+j]+q[n+j];
              }
            w[3*n+j]=w[j]+g*h;
            y[j]=w[3*n+j];
          }
        s=t+h;
        gtnr2f(s,y,n,d);
        for (j=0; j<=n-1; j++) q[3*n+j]=d[j];
        for (j=0; j<=n-1; j++)
          { if (p[j]<=0.0)
              { dy=q[j]+2.0*(q[n+j]+q[n+n+j]);
                dy=(dy+q[n+n+n+j])*h/6.0;
              }
            else
              { dy=-3.0*(q[j]+p[j]*w[j])+2.0*(q[n+j]
                          +p[j]*w[n+j]);
                dy=dy+2.0*(q[n+n+j]+p[j]*w[n+n+j]);
                dy=dy-(q[n+n+n+j]+p[j]*w[n+n+n+j]);
                dy=dy*r[n+n+j]+q[j]*r[n+j];
                dy1=q[j]-q[n+j]-q[n+n+j]+q[n+n+n+j];
                dy1=dy1+(w[j]-w[n+j]-w[n+n+j]+w[n+n+n+j])*p[j];
                dy=(dy+4.0*dy1*r[n+n+n+j])*h;
              }
            y[j]=w[j]+dy;
            z[j*k+i]=y[j];
          }
      }
    free(d); free(p); free(w); free(q); free(r);
    return;
}
/////////////////////////////////////////////////////////////
int ggear(double a,double b,double hmin,double hmax,double h,double eps,
		int n,double y0[],int k,double t[],double z[],
		void (*ggearf)(double t,double y[],int n,double d[]),
		void (*ggears)(double t,double y[],int n,double p[]))
{ 
	//extern void ggearf();
    //extern void ggears();
    //extern int brinv();
    int kf,jt,nn,nq,i,m,irt1,irt,j,nqd,idb;
    int iw,j1,j2,nt,nqw,l,*iis,*jjs;
    double aa[7],hw,hd,rm,t0,td,r,dd,pr1,pr2,pr3,rr;
    double enq1,enq2,enq3,eup,e,edwn,bnd,r1;
    static double pp[7][3]={ {2.0,3.0,1.0},{4.5,6.0,1.0},
           {7.333,9.167,0.5},{10.42,12.5,0.1667},
           {13.7,15.98,0.04133},{17.15,1.0,0.008267},
           {1.0,1.0,1.0}};
    double *d,*p,*s,*s02,*ym,*er,*yy,*y;
    d=(double*)malloc(n*sizeof(double));
    p=(double*)malloc(n*n*sizeof(double));
    s=(double*)malloc(10*n*sizeof(double));
    s02=(double*)malloc(n*sizeof(double));
    ym=(double*)malloc(n*sizeof(double));
    er=(double*)malloc(n*sizeof(double));
    yy=(double*)malloc(n*sizeof(double));
    y=(double*)malloc(8*n*sizeof(double));
    iis=(int*)malloc(n*sizeof(int));
    jjs=(int*)malloc(n*sizeof(int));
    aa[1]=-1.0; jt=0; nn=0; nq=1; t0=a;
    for (i=0; i<=8*n-1; i++) y[i]=0.0;
    for (i=0; i<=n-1; i++) 
       { y[i*8]=y0[i]; yy[i]=y[i*8];}
    ggearf(t0,yy,n,d);
    for (i=0; i<=n-1; i++) y[i*8+1]=h*d[i];
    hw=h; m=2;
    for (i=0; i<=n-1; i++) ym[i]=1.0;
    l20:
    irt=1; kf=1; nn=nn+1;
    t[nn-1]=t0;
    for (i=0; i<=n-1; i++) z[i*k+nn-1]=y[i*8];
    if ((t0>=b)||(nn==k))
      { free(d); free(p); free(s); free(s02);
        free(ym); free(er); free(yy); free(iis);
        free(jjs); free(y); return(kf);}
    for (i=0; i<=n-1; i++)
      for (j=0; j<=m-1; j++) s[i*10+j]=y[i*8+j];
    hd=hw;
    if (h!=hd)
      { rm=h/hd; irt1=0;
        rr=fabs(hmin/hd);
        if (rm<rr) rm=rr;
        rr=fabs(hmax/hd);
        if (rm>rr) rm=rr;
	r=1.0; irt1=irt1+1;
        for (j=1; j<=m-1; j++)
          { r=r*rm;
            for (i=0; i<=n-1; i++)
              y[i*8+j]=s[i*10+j]*r;
          }
        h=hd*rm;
        for (i=0; i<=n-1; i++)
          y[i*8]=s[i*10];
        idb=m;
      }
    nqd=nq; td=t0; rm=1.0;
    if (jt>0) goto l80;
    l60:
    switch (nq)
      { case 1: aa[0]=-1.0; break;
        case 2: aa[0]=-2.0/3.0; aa[2]=-1.0/3.0; break;
        case 3: aa[0]=-6.0/11.0; aa[2]=aa[0];
                aa[3]=-1.0/11.0; break;
        case 4: aa[0]=-0.48; aa[2]=-0.7; aa[3]=-0.2;
                aa[4]=-0.02; break;
        case 5: aa[0]=-120.0/274.0; aa[2]=-225.0/274.0;
                aa[3]=-85.0/274.0; aa[4]=-15.0/274.0;
                aa[5]=-1.0/274.0; break;
        case 6: aa[0]=-720.0/1764.0; aa[2]=-1624.0/1764.0;
                aa[3]=-735.0/1764.0; aa[4]=-175.0/1764.0;
                aa[5]=-21.0/1764.0; aa[6]=-1.0/1764.0;
                break;
        default: { free(d); free(p); free(s); free(s02);
                   free(ym); free(er); free(yy);
                   free(iis); free(jjs); free(y); return(-2);
                 }
      }
    m=nq+1; idb=m;
    enq2=0.5/(nq+1.0); enq3=0.5/(nq+2.0);
    enq1=0.5/(nq+0.0);
    eup=pp[nq-1][1]*eps; eup=eup*eup;
    e=pp[nq-1][0]*eps; e=e*e;
    edwn=pp[nq-1][2]*eps; edwn=edwn*edwn;
    if (edwn==0.0)
      { for (i=0; i<=n-1; i++)
          for (j=0; j<=m-1; j++)
            y[i*8+j]=s[i*10+j];
        h=hd; nq=nqd; jt=nq;
        free(d); free(p); free(s); free(s02);
        free(ym); free(er); free(yy); free(iis);
        free(jjs); free(y); return(-4);
      }
    bnd=eps*enq3/(n+0.0);
    iw=1;
    if (irt==2)
      { r1=1.0;
        for (j=1; j<=m-1; j++)
          { r1=r1*r;
            for (i=0; i<=n-1; i++)
              y[i*8+j]=y[i*8+j]*r1;
          }
        idb=m;
        for (i=0; i<=n-1; i++)
          if (ym[i]<fabs(y[i*8]))
            ym[i]=fabs(y[i*8]);
        jt=nq;
        goto l20;
      }
    l80:
    t0=t0+h;
    for (j=2; j<=m; j++)
      for (j1=j; j1<=m; j1++)
        { j2=m-j1+j-1;
          for (i=0; i<=n-1; i++)
            y[i*8+j2-1]=y[i*8+j2-1]+y[i*8+j2];
        }
    for (i=0; i<=n-1; i++) er[i]=0.0;
    j1=1; nt=1;
    for (l=0; l<=2; l++)
      { if ((j1!=0)&&(nt!=0))
          { for (i=0; i<=n-1; i++) yy[i]=y[i*8];
            ggearf(t0,yy,n,d);
            if (iw>=1)
              { for (i=0; i<=n-1; i++) yy[i]=y[i*8];
                ggears(t0,yy,n,p);
                r=aa[0]*h;
                for (i=0; i<=n-1; i++)
                  for (j=0; j<=n-1; j++)
                    p[i*n+j]=p[i*n+j]*r;
                for (i=0; i<=n-1; i++)
                  p[i*n+i]=1.0+p[i*n+i];
                iw=-1;
                jjs[0]=brinv(p,n);//#include "ch2.h"已經在ch1中存在
                j1=jjs[0];
              }
            if (jjs[0]!=0)
              { for (i=0; i<=n-1; i++)
                  s02[i]=y[i*8+1]-d[i]*h;