"$d_Qndfs""ODE runge-kutta 4, C code.txt"



#include "nrutil.h"

void rk4(y,dydx,n,x,h,yout,derivs)
float dydx[],h,x,y[],yout[];
int n;
void ( * derivs)();
{
	int i;
	float xh,hh,h6, * dym, * dyt, * yt;

	dym := vector(1,n);
	dyt := vector(1,n);
	yt := vector(1,n);
	hh := h * 0.5;
	h6 := h / 6.0;
	xh := x + hh;
	for (i := 1;i<=n;i +  + ) yt[i] := y[i] + hh * dydx[i];
	( * derivs)(xh,yt,dyt);
	for (i := 1;i<=n;i +  + ) yt[i] := y[i] + hh * dyt[i];
	( * derivs)(xh,yt,dym);
	for (i := 1;i<=n;i +  + ) {
		yt[i] := y[i] + h * dym[i];

	}
	( * derivs)(x + h,yt,dyt);
	for (i := 1;i<=n;i +  + )
		yout[i] := y[i] + h6 * (dydx[i] + dyt[i] + 2.0 * dym[i]);
	free_vector(yt,1,n);
	free_vector(dyt,1,n);
	free_vector(dym,1,n);
}



#include "nrutil.h"

float  *  * y, * xx;

void rkdumb(vstart,nvar,x1,x2,nstep,derivs)
float vstart[],x1,x2;
int nstep,nvar;
void ( * derivs)();
{
	void rk4();
	int i,k;
	float x,h;
	float  * v, * vout, * dv;

	v := vector(1,nvar);
	vout := vector(1,nvar);
	dv := vector(1,nvar);
	for (i := 1;i<=nvar;i +  + ) {
		v[i] := vstart[i];
		y[i][1] := v[i];
	}
	xx[1] := x1;
	x := x1;

	for (k := 1;k<=nstep;k +  + ) {
		( * derivs)(x,v,dv);
		rk4(v,dv,nvar,x,h,vout,derivs);
		if ((float)(x + h)  =  x) nrerror("Step size too small in routine rkdumb");

		xx[k + 1] := x;
		for (i := 1;i<=nvar;i +  + ) {
			v[i] := vout[i];
			y[i][k + 1] := v[i];
		}
	}
	free_vector(dv,1,nvar);
	free_vector(vout,1,nvar);
	free_vector(v,1,nvar);
}