#include <stdio.h>
#include <string.h>

#define A 0
#define B 10
#define K 1

double expected_u (double x, double t);
double initial_condition (double x); //t=0
double left_border_condition (double t); //x=A
double right_border_condition (double t); //x=B
double right_part_func (double x, double t);
void explicit_next_layer_calculate (double *prev_layer, double *next_layer, double layer_t, double x_step, double t_step);
void pr_ar (double * arr, int dim);

int main (void) {
  int input_param = -1;
  double end_time;
  double x_step, t_step;
  int x_dim, t_dim;
  while (input_param) {
    printf ("Please, type 1 if you want to use steps or type 2 if you want to use count of segments. Type 0 if you want to exit. >>");
    scanf ("%d", &input_param);
    if (input_param) {
      printf ("Type end_time. >>");
      scanf ("%lf", &end_time);
      switch (input_param)
      {
	case 1: {
	  printf ("Type range step (h). >>");
	  scanf ("%lf", &x_step);
	  printf ("Type time step (tau). >>");
	  scanf ("%lf", &t_step);
	  x_dim = (B - A) / x_step;
	  t_dim = end_time / t_step;
	  printf ("t_dim %d", t_dim);
	  if ((B - A - x_step * x_dim) > 0) {
	    x_dim += 2;
	  }
	  else {
	    x_dim += 1;
	  }
	  if ((end_time - t_step * t_dim)  > 0) {
	    t_dim += 2;
	  }
	  else {
	    t_dim += 1;
	  }
	  break;
	}
	case 2: {
	  printf ("Type count of segments for x axis. >>");
	  scanf ("%d", &x_dim);
	  printf ("Type count of segments for t axis. >>");
	  scanf ("%d", &t_dim);
	  x_step = (B - A) / x_dim;
	  t_step = end_time / t_dim;
	  x_dim ++;
	  t_dim ++;
	  break;
	}
	default: {
	  printf ("Invalid command. Type 1 or 2 to continue. Type 0 to exit. >>");
	  scanf ("%d", &input_param);
	}
      }
      printf ("Input completed. Following values were typed:\n\th = %.3lf\n\ttau = %.3lf\n", x_step, t_step);
      printf ("Count of points per each axis:\n\tfor x: %d\n\tfor t: %d\n", x_dim, t_dim);


      printf("x_dim %d", x_dim);
      double prev_layer [x_dim];
      double next_layer [x_dim];
      printf("wtf");
      for (int i = 0; i < x_dim; i++) {
	prev_layer[i] = initial_condition (A + x_step * i);
	//printf("%d\t",i);
      }

      double t = t_step;
      while (t < end_time) {
	explicit_next_layer_calculate (prev_layer, next_layer, t, x_step, t_step);
	memcpy (prev_layer, next_layer, x_dim*(sizeof(double)));
	t += t_step;
      }

      if (t-t_step > 0.0) {
	explicit_next_layer_calculate (prev_layer, next_layer, end_time, x_step, (B - t - t_step));
      }
      else {
	explicit_next_layer_calculate (prev_layer, next_layer, end_time, x_step, t_step);
      }

      pr_ar(next_layer, x_dim);
    }
  }

  return 0;
}

void explicit_next_layer_calculate (double *prev_layer, double *next_layer, double layer_t, double x_step, double t_step) {

  next_layer[0] = left_border_condition(layer_t);
  int i = 1;
  int x = A + t_step;
  while (x + x_step < B) {
    next_layer[i] = prev_layer[i-1] + (K * t_step / (x_step * x_step)) * (prev_layer[i-1] - 2 * prev_layer[i] + prev_layer[i+1]) +
      t_step * right_part_func(x,layer_t);
    i++;
    x = x + x_step;
  }
  next_layer[i] = right_border_condition(layer_t);
}

double expected_u (double x, double t) {
  return x * x * x + 1 + 4 * t * t;
}

double initial_condition (double x){ //t=0
  return x * x * x +1;
}

double left_border_condition (double t) { //x=A
  return A * A * A + 1 + 4 * t * t;
}

double right_border_condition (double t) { //x=B
  return B * B * B + 1 + 4 * t * t;
}

double right_part_func (double x, double t) {
  return 8 * t - 6 * x;
}

void pr_ar (double * arr, int dim) {
  for (int i = 0; i < dim; i++) {
    printf ("%.3lf/t", arr[i]);
  }
}