#include <math.h>#define NRANSI#include "nrutil.h"#define SAFETY 0.9#define GROW 1.5#define PGROW -0.25#define SHRNK 0.5#define PSHRNK (-1.0/3.0)#define ERRCON 0.1296#define MAXTRY 40#define GAM (1.0/2.0)#define A21 2.0#define A31 (48.0/25.0)#define A32 (6.0/25.0)#define C21 -8.0#define C31 (372.0/25.0)#define C32 (12.0/5.0)#define C41 (-112.0/125.0)#define C42 (-54.0/125.0)#define C43 (-2.0/5.0)#define B1 (19.0/9.0)#define B2 (1.0/2.0)#define B3 (25.0/108.0)#define B4 (125.0/108.0)#define E1 (17.0/54.0)#define E2 (7.0/36.0)#define E3 0.0#define E4 (125.0/108.0)#define C1X (1.0/2.0)#define C2X (-3.0/2.0)#define C3X (121.0/50.0)#define C4X (29.0/250.0)#define A2X 1.0#define A3X (3.0/5.0)void stiff(float y[], float dydx[], int n, float *x, float htry, float eps,	float yscal[], float *hdid, float *hnext,	void (*derivs)(float, float [], float [])){	void jacobn(float x, float y[], float dfdx[], float **dfdy, int n);	void lubksb(float **a, int n, int *indx, float b[]);	void ludcmp(float **a, int n, int *indx, float *d);	int i,j,jtry,*indx;	float d,errmax,h,xsav,**a,*dfdx,**dfdy,*dysav,*err;	float *g1,*g2,*g3,*g4,*ysav;	indx=ivector(1,n);	a=matrix(1,n,1,n);	dfdx=vector(1,n);	dfdy=matrix(1,n,1,n);	dysav=vector(1,n);	err=vector(1,n);	g1=vector(1,n);	g2=vector(1,n);	g3=vector(1,n);	g4=vector(1,n);	ysav=vector(1,n);	xsav=(*x);	for (i=1;i<=n;i++) {		ysav[i]=y[i];		dysav[i]=dydx[i];	}	jacobn(xsav,ysav,dfdx,dfdy,n);	h=htry;	for (jtry=1;jtry<=MAXTRY;jtry++) {		for (i=1;i<=n;i++) {			for (j=1;j<=n;j++) a[i][j] = -dfdy[i][j];			a[i][i] += 1.0/(GAM*h);		}		ludcmp(a,n,indx,&d);		for (i=1;i<=n;i++)			g1[i]=dysav[i]+h*C1X*dfdx[i];		lubksb(a,n,indx,g1);		for (i=1;i<=n;i++)			y[i]=ysav[i]+A21*g1[i];		*x=xsav+A2X*h;		(*derivs)(*x,y,dydx);		for (i=1;i<=n;i++)			g2[i]=dydx[i]+h*C2X*dfdx[i]+C21*g1[i]/h;		lubksb(a,n,indx,g2);		for (i=1;i<=n;i++)			y[i]=ysav[i]+A31*g1[i]+A32*g2[i];		*x=xsav+A3X*h;		(*derivs)(*x,y,dydx);		for (i=1;i<=n;i++)			g3[i]=dydx[i]+h*C3X*dfdx[i]+(C31*g1[i]+C32*g2[i])/h;		lubksb(a,n,indx,g3);		for (i=1;i<=n;i++)			g4[i]=dydx[i]+h*C4X*dfdx[i]+(C41*g1[i]+C42*g2[i]+C43*g3[i])/h;		lubksb(a,n,indx,g4);		for (i=1;i<=n;i++) {			y[i]=ysav[i]+B1*g1[i]+B2*g2[i]+B3*g3[i]+B4*g4[i];			err[i]=E1*g1[i]+E2*g2[i]+E3*g3[i]+E4*g4[i];		}		*x=xsav+h;		if (*x == xsav) nrerror("stepsize not significant in stiff");		errmax=0.0;		for (i=1;i<=n;i++) errmax=FMAX(errmax,fabs(err[i]/yscal[i]));		errmax /= eps;		if (errmax <= 1.0) {			*hdid=h;			*hnext=(errmax > ERRCON ? SAFETY*h*pow(errmax,PGROW) : GROW*h);			free_vector(ysav,1,n);			free_vector(g4,1,n);			free_vector(g3,1,n);			free_vector(g2,1,n);			free_vector(g1,1,n);			free_vector(err,1,n);			free_vector(dysav,1,n);			free_matrix(dfdy,1,n,1,n);			free_vector(dfdx,1,n);			free_matrix(a,1,n,1,n);			free_ivector(indx,1,n);			return;		} else {			*hnext=SAFETY*h*pow(errmax,PSHRNK);			h=(h >= 0.0 ? FMAX(*hnext,SHRNK*h) : FMIN(*hnext,SHRNK*h));		}	}	nrerror("exceeded MAXTRY in stiff");}#undef SAFETY#undef GROW#undef PGROW#undef SHRNK#undef PSHRNK#undef ERRCON#undef MAXTRY#undef GAM#undef A21#undef A31#undef A32#undef C21#undef C31#undef C32#undef C41#undef C42#undef C43#undef B1#undef B2#undef B3#undef B4#undef E1#undef E2#undef E3#undef E4#undef C1X#undef C2X#undef C3X#undef C4X#undef A2X#undef A3X#undef NRANSI