36 lines
654 B
C
36 lines
654 B
C
|
|
#include "nrutil.h"
|
|
|
|
void stoerm(y,d2y,nv,xs,htot,nstep,yout,derivs)
|
|
float d2y[],htot,xs,y[],yout[];
|
|
int nstep,nv;
|
|
void (*derivs)();
|
|
{
|
|
int i,n,neqns,nn;
|
|
float h,h2,halfh,x,*ytemp;
|
|
|
|
ytemp=vector(1,nv);
|
|
h=htot/nstep;
|
|
halfh=0.5*h;
|
|
neqns=nv/2;
|
|
for (i=1;i<=neqns;i++) {
|
|
n=neqns+i;
|
|
ytemp[i]=y[i]+(ytemp[n]=h*(y[n]+halfh*d2y[i]));
|
|
}
|
|
x=xs+h;
|
|
(*derivs)(x,ytemp,yout);
|
|
h2=h*h;
|
|
for (nn=2;nn<=nstep;nn++) {
|
|
for (i=1;i<=neqns;i++)
|
|
ytemp[i] += (ytemp[(n=neqns+i)] += h2*yout[i]);
|
|
x += h;
|
|
(*derivs)(x,ytemp,yout);
|
|
}
|
|
for (i=1;i<=neqns;i++) {
|
|
n=neqns+i;
|
|
yout[n]=ytemp[n]/h+halfh*yout[i];
|
|
yout[i]=ytemp[i];
|
|
}
|
|
free_vector(ytemp,1,nv);
|
|
}
|