// Crank-Niclson 1st Order ODE solver
#include <stdio.h>
#include <math.h>

double f(double t, double y);
double df(double t, double y);
double analyticalsol(double x);
void crankniclson(double dt, double t0, double y0, double tmax);
void plotfunction(double x0, double x1, int N);

// Driver program
int main() 
{ 
    // Initial values
    crankniclson(0.125,0,1,50); 
    plotfunction(0, 50,1000);
    return 0; 
} 


void crankniclson(double dt, double t0, double y0, double tmax) 
{ 
    int i;
   double y,yy,t;
   FILE *output;

   output = fopen("crankniclson.dat", "w");//opening file. 

    t=t0;
    y=y0;
    fprintf(output,"%.15f\t%.15f\n",t0,y0);

    while (t <= tmax) 
    {
       t=t+dt;
       yy=y;
       // Netwton Raphson iterations
       for(i=1;i<1000;i++)
       {y=y-(-y+yy+0.5*dt*f(t-dt,yy)+0.5*dt*f(t,y))/(-1+0.5*dt*df(t,y));
       }

       y=yy+dt*f(t,y); 
      fprintf(output,"%.15f\t%.15f\n",t,y);

    } 
     fclose(output);//closing file.

} 


void plotfunction(double x0, double x1,int N)
{ int i;
  double x,dx;
  FILE *output;

      dx=(x1-x0)/N;


output = fopen("function.dat", "w");//opening file. 
 for(i=0;i<N;i++){
    x=x0+dx*i;
fprintf(output,"%.15f\t%.15f\n",x,analyticalsol(x));
 }
 
 fclose(output);//closing file.

}


// Define function

double f(double t, double y) 
{ 
    double fun;

fun=-y;
    return fun;
} 

double df(double t, double y) 
{ 
    double fun;

fun=-1;
    return fun;
} 
double analyticalsol(double x) 
{ 
    double y;

y=exp(-x);
    return y;
}