diff --git a/Makefile b/Makefile
index 6dde9cf72e8a8006f2f80b1885d1c923dc98b3bd..7bdc4a4672c81648abeb04a7709273a332431c4f 100644
--- a/Makefile
+++ b/Makefile
@@ -7,7 +7,7 @@ PROJNAME = miscmaths
USRINCFLAGS = -I${INC_NEWMAT} -I${INC_PROB} -I${INC_ZLIB}
USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_PROB} -L${LIB_ZLIB}
-OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o sparse_matrix.o sparsefn.o
+OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o sparse_matrix.o sparsefn.o rungekutta.o
#OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o
LIBS = -lutils -lnewmat -lprob -lm
diff --git a/rungekutta.cc b/rungekutta.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8a14621bec93f6b6805b21dea6e0dff9cda3d5c2
--- /dev/null
+++ b/rungekutta.cc
@@ -0,0 +1,167 @@
+/* rungekutta.cc
+
+ Mark Woolrich, FMRIB Image Analysis Group
+
+ Copyright (C) 2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "rungekutta.h"
+
+namespace MISCMATHS {
+
+
+void rk(ColumnVector& ret, const ColumnVector& y, const ColumnVector& dy, float x, float h, const Derivative& deriv,const ColumnVector& paramvalues)
+{
+ Tracer tr("rk");
+
+ float hh=h*0.5;
+ float xh=x+hh;
+
+ //first step
+ ColumnVector yt=y+hh*dy;
+
+ //second step
+ ColumnVector dyt = deriv.evaluate(xh,yt,paramvalues);
+ yt=y+hh*dyt;
+
+ //third step
+ ColumnVector dym = deriv.evaluate(xh,yt,paramvalues);
+ yt=y+h*dym;
+ dym=dym+dyt;
+
+ //fourth step
+ dyt = deriv.evaluate(x+h,yt,paramvalues);
+
+ //addup
+ ret = y+h*(dy+dyt+2*dym)/6;
+}
+
+void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry, float eps, const Derivative& deriv,const ColumnVector& paramvalues)
+{
+ Tracer tr("rkqc");
+
+ float xsav = x;
+ ColumnVector dysav = dy;
+ ColumnVector ysav = y;
+ float h = htry;
+ float hdid;
+ ColumnVector ytemp;
+
+ while(true)
+ {
+ // take 2 1/2 step sizes
+
+ // first 1/2 step
+ float hh=h*0.5;
+
+ rk(ytemp,ysav,dysav,xsav,hh,deriv,paramvalues);
+
+ // second 1/2 step
+ x=xsav+hh;
+ dy = deriv.evaluate(x,ytemp,paramvalues);
+ rk(y,ytemp,dysav,xsav,hh,deriv,paramvalues);
+ x=xsav+h;
+ if(x==xsav) cerr << "step size too small" << endl;
+
+ // take large step size
+ rk(ytemp,ysav,dysav,xsav,h,deriv,paramvalues);
+
+ // eval accuracy
+ float errmax = 0.0;
+ for(int i=1; i<=y.Nrows(); i++)
+ {
+ //errmax=max(abs((y-ytemp)./y));
+
+ float tmp = fabs((y(i)-ytemp(i))/y(i));
+ if(tmp > errmax) errmax = tmp;
+ }
+
+ errmax=errmax/eps;
+
+ if(errmax <=1.0)
+ {
+ // step OK, compute step size for next step
+ hdid=h;
+
+ if(errmax>6e-4)
+ hnext=h*std::exp(-0.2*std::log(errmax));
+ else
+ hnext=4*h;
+
+ break;
+ }
+ else
+ {
+ // step too large,
+ h=h*std::exp(-0.25*std::log(errmax));
+ }
+ }
+
+ y = y+(y-ytemp)/15;
+}
+
+void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVector& ystart, float x1, float x2, float eps, float hmin, const Derivative& deriv,const ColumnVector& paramvalues)
+{
+ Tracer tr("runge_kutta");
+
+ int MAXSTEP=1000;
+
+ ColumnVector y = ystart;
+
+ float x=x1;
+ xp.ReSize(MAXSTEP,1);
+ xp = 0;
+ xp(1) =x1;
+
+ float h=hp(1);
+ hp.ReSize(MAXSTEP,1);
+ hp = 0;
+
+ yp.ReSize(MAXSTEP,y.Nrows());
+ yp = 0;
+
+ int kout=1;
+
+ ColumnVector dy;
+
+ for(int k=1; k <= MAXSTEP; k++)
+ {
+ dy = deriv.evaluate(x,y,paramvalues);
+
+ // store results:
+ xp(kout)=x;
+ yp.Row(kout)=y;
+ hp(kout)=h;
+
+ kout=kout+1;
+
+ // stop overshoot of step past x2:
+ if((x+h-x2)*(x+h-x1)>0) h=x2-x;
+
+ float hnext = 0.0;
+ rkqc(y,x,hnext,dy,h,eps,deriv,paramvalues);
+
+ if((x-x2)*(x2-x1) >= 0.0)
+ {
+ xp(kout)=x;
+ yp.Row(kout)=y;
+ hp(kout)=h;
+ //kout=kout+1;
+
+ xp = xp.Rows(1,kout);
+ yp = yp.Rows(1,kout);
+
+ return;
+ }
+ else
+ {
+ if(hnext<=hmin) cerr << "step size too small" << endl;
+ h=hnext;
+ }
+
+ }
+ cerr << "too many steps" << endl;
+}
+
+}
diff --git a/rungekutta.h b/rungekutta.h
new file mode 100644
index 0000000000000000000000000000000000000000..2f57123e8dabbd19029781768185ce879aff1747
--- /dev/null
+++ b/rungekutta.h
@@ -0,0 +1,48 @@
+/* rungekutta.h
+
+ Mark Woolrich - FMRIB Image Analysis Group
+
+ Copyright (C) 2002 University of Oxford */
+
+/* COPYRIGHT */
+
+#if !defined(rungekutta_h)
+#define rungekutta_h
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include "newimage/newimageall.h"
+
+using namespace NEWIMAGE;
+
+namespace MISCMATHS {
+
+class Derivative
+{
+public:
+ Derivative(int pny) : ny(pny), dy(pny) {}
+
+ // x is time point to evaluate at
+ // y is state variables
+ // paramvalues are "constants" in the diff eqn
+ virtual const ColumnVector& evaluate(float x,const ColumnVector& y,const ColumnVector& paramvalues) const = 0;
+
+ virtual ~Derivative(){};
+
+protected:
+ int ny;
+ mutable ColumnVector dy;
+};
+
+void rk(ColumnVector& ret, const ColumnVector& y, const ColumnVector& dy, float x, float h, const Derivative& deriv,const ColumnVector& paramvalues);
+
+void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry, float eps, const Derivative& deriv,const ColumnVector& paramvalues);
+
+void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVector& ystart, float x1, float x2, float eps, float hmin, const Derivative& deriv,const ColumnVector& paramvalues);
+
+}
+
+#endif