//cc dumka3.c dumkaC.c -lm
#include <stdio.h>
#include "math.h"
FILE *solution;

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*		FitzHugh and Nagumo problem
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*
*   This program solves a system of ODEs resulting from the 
*    space discretization of 
*   the FitzHugh and Nagumo equations (u=u(x,t),v=v(x,t)):
*    
*   u_t=u_{xx}-f(u)-v  for 0<= t <=400, 0<= x <= 100
*   v_t=n(u-b*v)       where a=0.139, n=0.008, b=2.54
*
*    with initial condition 
*
*    u(x,0)=v(x,0)=0
*
*   and boundary conditions
*
*    u_t(0,t)=-0.3,  u_t(100,t)=0
*
*    The function f is defined by 
*
*    f(u)=u*(u-a)*(u-1)
*
*    We discretize the space variables with
*    x_i=(2*i+1)/4, for i=0,1,...,199.
*    We obtain a system of 400 equations.
*    The spectral radius of the Jacobian can
*    be estimated with the Gerhgorin theorem   
*    Thus we provide an external function COUR,
*    giving the spectral radius of the Jacobian matrix.
*    As output point we choose t_out=400 
*
* The algorithm dumka3 and dumka4 is described in the article:
* Medovikov A.A. High order explicit methods for parabolic 
*                equations. BIT, V38,No2,pp.372-390,1998 
*
* Many thanks Prof. Lebedev V.I., Prof. Wanner G. and Prof. Hairer E.
* and Prof. A. Abdulle
*/

//parameters of the problem
double alpha = 0.139e0;
double beta = 2.54e0;
double eta = 0.008e0;
double rhop = 0.3e0;
double dif = 0.010e0;
double tau, rhodn, dis;

void rhs(long neqn, double* y, double t, double* result ) {
	int i;
	// ---------- discretisation  -------
	double fv=((y[0]-alpha-1.0e0)*y[0]+alpha)*y[0];
	result[0]=dis*(rhodn-y[0]+y[2])-fv-y[1];
	for ( i=2; i < neqn-3; i=i+2 ) {
		fv=((y[i]-alpha-1.0e0)*y[i]+alpha)*y[i];
		result[i]=dis*(y[i-2]-2.0e0*y[i]+y[i+2])-fv-y[i+1];
	}
	fv=((y[neqn-2]-alpha-1.0e0)*y[neqn-2]+alpha)*y[neqn-2];
	result[neqn-2]=dis*(y[neqn-4]-y[neqn-2])-fv-y[neqn-1];
	// ----
	for ( i=1; i < neqn; i=i+2 ) {
		result[i]=eta*y[i-1]+tau*y[i];
	}
};
/*
*    The subroutine cour gives an estimation of the spectral 
*    radius of the Jacobian matrix of the problem rho and returns
*    2/rho, which is tha maximum possible step in explicit Euler method.
-
*/			
double cour(double* y, double t, 
			void (f(long neqn, double* y, double t, double* result))) {
				double rho=2.e0/17.e0;
				return rho;
			};
int main(int argc, char** argv)
{
	void dumka3(long neqn, double* t, const double tend, const double h0,
		const double atol, const double rtol, 
		void (f(long neqn, double* y, double t, double* result)),
		double (cour(double* y, double t, void (f(long neqn, double* y, double t, double* result)))),
		double* y,double* z0,double* z1,double* z2,
		double* z3,double* oldEigenVector,
		int isComputeEigenValue);
	double y[400],z0[400], z1 [400],z2[400],z3[400], oldEigenVector[400];
	int isComputeEigenValue =(1==0);
	long neqn = 400;
	int i;
	double  t = 0.e0, tend = 400.e0, h0 = 0.0001e0, atol = 0.01e0, rtol = 0.01e0;
	double hd=200e0;
    dis=(hd*dif)*(hd*dif);
    tau=-eta*beta;
    rhodn=rhop/(hd*dif);

	// ----- initial values -----
	for ( i=0; i < neqn; i++ ) {
		y[i]=0.e-0;
	}; 

	dumka3(neqn, &t,  tend, h0,
		atol, rtol, rhs, cour, 
		y,z0,z1,z2, z3, oldEigenVector, isComputeEigenValue);

	// output result
	solution = fopen( "u", "w+" );
	for ( i=0; i < 400; i=i+2 ) {
        fprintf( solution, "%e \n",  y[i] );
	}
	fclose( solution );
	solution = fopen( "v", "w+" );
	for ( i=1; i < 400; i=i+2 ) {
        fprintf( solution, "%e \n",  y[i] );
	}
	fclose( solution );

	return 0;
};