#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct {
	double 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;
    //FILE *fp2=fopen("cosh.csv","w");
	for (int j = 0; j <= jmax; j++) {

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

	}

}
void ujeins(soli_param *sp) {
	double jmax = sp->jmax;
	double d = sp->d;
	double h = sp->h;
    double x=sp->xmin;
	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*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*u0[j+1]) / (2 * pow(h, 3))) + u0[j];

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

		} else if (j == jmax-1) {
			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]+ 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];
    double t=sp->tmin;
    double x=sp->xmin;
    int n;
    
	for ( n = 2; 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];*/
                unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(un[j+2]-2*un[j+1]+2*un[j-1]) ) +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];*/
                unp[j]=(2*(d/h) *(un[j+1]+un[j])*(un[j+1]))-(d/(h*h*h)*(un[j+2]-2*un[j+1]) ) +unm[j];  

			} else if (j == jmax - 2) {
				/*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];*/
                unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(-2*un[j+1]+2*un[j-1]-un[j-2]) ) +unm[j];  

			} else if (j == jmax - 1) {
				/*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];*/
                unp[j]=(2*(d/h) *(un[j]+un[j-1])*(-un[j-1]))-(d/((h*h*h))*(2*un[j-1]-un[j-2]) ) +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];*/
                                                
                 unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(un[j+2]-2*un[j+1]+2*un[j-1]-un[j-2]) ) +unm[j];                     
			}
			//unm[j]=un[j];
            //un[j]=unp[j];
			

		}
		for (int j = 0; j < jmax; j++) {
    		  unm[j]=un[j];
    		  un[j]=unp[j];
            }
		
		
        
        

	}
	FILE *fp= fopen("numplot8.csv", "w");
	printf("%f\t%f\n", d*n, d*nmax);
	for(int j=0;j<jmax;j++){
		
		fprintf(fp, "%f\t%f\n",x,un[j]);
		//printf( "%f\t%f\n",x,un[j]);
        x=x+h;

	}
	fclose(fp);
}


int main(void) {
    soli_param sp;
    sp.N = 1;
    sp.h=0.1;
    sp.d=0.0001;
    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;
    //sp.nmax=500;
    printf("%d\n", sp.nmax);
    u0 = (double*) malloc(sp.jmax * sizeof(double));
    ustart(&sp);
    uj1 = (double*) malloc(sp.jmax * sizeof(double));
    ujeins(&sp);
    nextu(&sp);
    free(u0);
    free(uj1);

}