#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct {
	int N;
	double d;
	double h;
	double xmin, xmax, tmin, tmax;
	double j;
	double n;
	int jmax;
	int nmax;

} soli_param;
double *u0;
double *uj1;
void ustart(soli_param *sp) {
	double N = sp->N;
	double jmax = sp->jmax;
	double h = sp->h;
	double x = sp->xmin;

	for (int j = 0; j <= jmax; j++) {

		u0[j] = -N * (N + 1) / (pow(cosh(x), 2));
		x = x + h;
		printf("%f\n", u0[j]);

	}

}
void ujeins(soli_param *sp) {
	double jmax = sp->jmax;
	double d = sp->n;
	double h = sp->h;
	for (int j = 0; j <= jmax; j++) {
		if (j == 1) {
			uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
					- (d * (u0[j + 2] + 2 * u0[j - 1]) / (2 * pow(h, 3)))
					+ u0[j];

		} else if (j == 0) {
			uj1[j] = (d * 6 * u0[j] * (u0[j + 1]) / (2 * h))
					- (d * (u0[j + 2]) / (2 * pow(h, 3))) + u0[j];

		} else if (j == jmax - 1) {
			uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
					- (d * (2 * u0[j - 1] - u0[j - 2]) / (2 * pow(h, 3)))
					+ u0[j];

		} else if (j == jmax) {
			uj1[j] = (d * 6 * u0[j] * (-u0[j - 1]) / (2 * h))
					- (d * (2 * u0[j - 1] - u0[j - 2]) / (2 * pow(h, 3)))
					+ u0[j];
		} else {
			uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
					- (d * (u0[j + 2] + 2 * u0[j - 1] - u0[j - 2])
							/ (2 * pow(h, 3))) + u0[j];
		}

	}

}
void nextu(soli_param *sp) {
	int jmax = sp->jmax;
	int nmax = sp->nmax;
	double d = sp->d;
	double h = sp->h;
	double * unm = u0;
	double * un = uj1;
	double unp[jmax];
	for (int n = 0; n <= nmax; n++) {
		for (int j = 0; j <= jmax; j++) {
			if (j == 1) {
				unp[j] = 2*d*(
								(un[j + 1] + un[j] + un[j - 1])
										* (un[j + 1] - un[j - 1]) / h
										- (un[j + 2] - 2 * un[j + 1]
												+ 2 * un[j - 1])
												/ (2 * pow(h, 3))) + unm[j];

			} else if (j == 0) {
				unp[j] = 2
						* d*(
								(un[j + 1] + un[j]) * (un[j + 1]) / h
										- (un[j + 2] - 2 * un[j + 1])
												/ (2 * pow(h, 3))) + unm[j];

			} else if (j == jmax - 1) {
				unp[j] = 2
						* d*(
								(un[j + 1] + un[j] + un[j - 1])
										* (un[j + 1] - un[j - 1]) / h
										- (-2 * un[j + 1] + 2 * un[j - 1]
												- un[j - 2]) / (2 * pow(h, 3)))
						+ unm[j];

			} else if (j == jmax) {
				unp[j] = 2
						* d*(
								(un[j] + un[j - 1]) * (-un[j - 1]) / h
										- (2 * un[j - 1] - un[j - 2])
												/ (2 * pow(h, 3))) + unm[j];

			} else {
				unp[j] = 2
						* d*(
								(un[j + 1] + un[j] + un[j - 1])
										* (un[j + 1] - un[j - 1]) / h
										- (un[j + 2] - 2 * un[j + 1]
												+ 2 * un[j - 1] - un[j - 2])
												/ (2 * pow(h, 3))) + unm[j];
			}


		}
		unm=un;
		un=unp;

	}
	FILE *fp= fopen("numplot2.csv", "w");
	double x=sp->xmin;
	for(int j=0;j<=jmax;j++){
		x=x+h;
		fprintf(fp, "%f\t%f\n",x,un[j]);
		printf( "%f\t%f\n",x,un[j]);

	}
	fclose(fp);
}


int main(void) {
soli_param sp;
sp.N = 2;
sp.h=0.1;
sp.d=pow(sp.h,3)/3;
sp.tmax=1;
sp.tmin=0;
sp.xmin=-30;
sp.xmax=30;
sp.jmax=(sp.xmax-sp.xmin)/sp.h;
sp.nmax=(sp.tmax-sp.tmin)/sp.d;
u0 = (double*) malloc(sp.jmax * sizeof(double));
ustart(&sp);
uj1 = (double*) malloc(sp.jmax * sizeof(double));
ujeins(&sp);
nextu(&sp);

}