Solution for 1 D heat conduction with Finite Difference not reducing error with increase in grid points

Question

I am solving 1 D heat conduction equation.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void grid(int nx, double xst, double xen, double *x, double *dx)
{
  int i;

  *dx = (xen-xst)/(double)(nx-1);

  for(i=0; i<nx; i++)
    x[i] = (double)i * (*dx); // ensure x[0] == 0.0 and x[nx-1] == 1.0
}

void get_exact_solution(int nx, double time, double *x, double *Texact)
{
// calculates exact solution
  int i;

  for(i=0; i<nx; i++)
    Texact[i] = erf((x[i]-0.5)/2.0/sqrt(time));
}


void enforce_bcs(int nx, double *x, double *T)
{
  T[0] = -1.0;
  T[nx-1] = 1.0;
}

void set_initial_condition(int nx, double *x, double *T)
{
  int i;

  for(i=0; i<nx; i++)
    T[i] = tanh((x[i]-0.5)/0.1);

  enforce_bcs(nx,x,T); //ensure BCs are satisfied at t = 0
}

void get_rhs(int nx, double dx, double *x, double *T, double *rhs)
{
  int i;
  double dxsq = dx*dx;

  // compute rhs. For this problem, d2T/dx2
  for(i=1; i<nx-1; i++)
    rhs[i] = (T[i+1]+T[i-1]-2.0*T[i])/dxsq;
}

void timestep_Euler(int nx, double dt, double dx, double *x, double *T, double *rhs)
{

  int i;

  // compute rhs
  get_rhs(nx,dx,x,T,rhs);

  // (Forward) Euler scheme
  for(i=1; i<nx-1; i++)
    T[i] = T[i] + dt*rhs[i];  // T^(it+1)[i] = T^(it)[i] + dt * rhs[i];

  // set Dirichlet BCs
  enforce_bcs(nx,x,T);

}

void output_soln(int nx, int it, double tcurr, double *x, double *T)
{
  int i;
  FILE* fp;
  char fname[100];

  sprintf(fname, "T_x_%04d.dat", it);
  //printf("\n%s\n", fname);

  fp = fopen(fname, "w");
  for(i=0; i<nx; i++)
    fprintf(fp, "%lf %lf\n", x[i], T[i]);
  fclose(fp);
}

int main()
{

  int nx;
  double *x, *T, *rhs, tst, ten, xst, xen, dx, dt, tcurr, *Texact;
  int i, it, num_time_steps, it_print;
  FILE* fp;  


  nx =10; //grid points

  xst = 0, //starting and ending in space
  xen = 1;
  tst=0; // start and end in time
  ten=1


  printf("Inputs are: %d %lf %lf %lf %lf\n", nx, xst, xen, tst, ten);

  x = (double *)malloc(nx*sizeof(double));
  T = (double *)malloc(nx*sizeof(double));
  Texact = (double *)malloc(nx*sizeof(double));
  rhs = (double *)malloc(nx*sizeof(double));


  grid(nx,xst,xen,x,&dx);         // initialize the grid

  set_initial_condition(nx,x,T);  // initial condition

  // prepare for time loop
  dt = 0.001; 

  num_time_steps = (int)((ten-tst)/dt) + 1; // why add 1 to this?
  it_print = num_time_steps/10;             // write out approximately 10 intermediate results

  // start time stepping loop
  for(it=0; it<num_time_steps; it++)
  {
    tcurr = tst + (double)it * dt;

    timestep_Euler(nx,dt,dx,x,T,rhs);    // update T

    // output soln every it_print time steps
    //if(it%it_print==0)
      //output_soln(nx,it,tcurr,x,T);
  }

  get_exact_solution( nx, ten, x, Texact);

  int loop;
  double sq_diff,total,err;
  total = 0;

  for(loop = 0; loop < nx; loop++)
      printf("%f ", T[loop]);

  printf("\n");

  for(loop = 0; loop < nx; loop++)
      printf("%f ", Texact[loop]);

  for(loop = 0; loop < nx; loop++){
      sq_diff = (T[loop] - Texact[loop])*(T[loop] - Texact[loop]);
      total = total+sq_diff;
    }

  err = sqrt(total);
  printf("\n L2 Norm error = %f",err);


  printf("\n");
  free(rhs);
  free(T);
  free(x);

  return 0;
}

Ideally, increasing number of grid point (nx in the program) should decrease the Error that I have calculated on final Time ten using exact solution Texact and numerical solution T.

The error should decrease with small increase in nx till the solution is stable and follows stability criteria (Von Neumann stability criteria.). But in my case this is not happening and error is increasing with increase in nx.

I am calculating the L2_norm error between T and Texact. Texact is calculated by the function get_exact_solution.

This plot shows the problem for different cases.

Please help in correcting this program.