Code:

double f(double Axt, double Mxt, double r, double k2, double u2) /* Definition of A(x,t) with M(x,t) carry in as parameter*/
{
double function;
function = (r*(1-(Axt/k2))*Mxt)-(u2+gamma)*Axt;
return(function);
}
double m(double Axt, double Mxt, double m1,double m2, double m3, double t, double x, double D, double v, double k1, double gamma, double u1) /* Definition of M(x,t) after explicit finite difference with A(x,t) carry in as parameter*/
{
double function1,function2,function3,function4,final;
function1 = (((t*D)/(x*x))-(t/x)*v*Mxt)*m1;
function2 = (1-((2*t*D)/(x*x)))+(((t/x)*(v*Mxt))*m2);
function3 = ((t*D)/(x*x))*m3;
function4 = ((gamma*t*Axt)*(1-(Mxt/k1)))-(t*u1*Mxt);
final = function1+function2+function3+function4;
return(final);
}
double Axt_runged (double Axt, double Mxt) // runge kutta A(x,t)
{
int N=100;
int i=0;
double t, y;
double a=0,b=1; // Initial values
double h=(b-a)/N; // Step size
double k1,k2,k3,k4;
t=a;
y=1;
while (i <= N)
{
k1=h*f(t,y, r,k2,u2);
k2=h*f(t+h/2,y+k1/2,r,k2,u2);
k3=h*f(t+h/2,y+k2/2,r,k2,u2);
k4=h*f(t+h,y+k3,r,k2,u2);
y=y+(k1+2*k2+2*k3+k4)/6; // Linear approximation to
i++; // Increments counter k by 1
t=t+h; // Increments t by the step size
}
getch();
return(y);
}