//Using shooting method to solve boundary value problem of the form y''=f(y',y,y); y(t0)=y0 and y(t1)=y1

// Shooting method for solving boundary value problem of the form d^2y/dx^2=f(dy/dx,y,t)

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

double f(double t, double y, double dy);
double analyticalsol(double x);
void shootingmethod(int N, double t0, double t1, double y0, double y1);
void plotfunction(double x0, double x1, int N);

// Driver program
int main() 
{ 
    // solve the problem defined by the function f
    // shootingmehtod(N,x0,x1,y0,y1)
    shootingmethod(1000,1.0,2.0,0.0,0.693); 
    plotfunction(1.0, 2.0,500);
    return 0; 
} 


// shooting method
void shootingmethod(int N, double t0, double t1, double y0, double y1) 
{ 
   double dalpha,alpha,y,u,k1,k2,k3,k4,t;
   double l1,l2,l3,l4,error,dt,yy1,yy2,alpha1,alpha2;
   int i;
   FILE *output;

    output = fopen("shootingmethod.dat", "w");//opening file. 
    dt=(t1-t0)/N;
    alpha=0.3;
    dalpha=0.3;
    t=t0;
    y=y0;
    u=alpha;
    
    fprintf(output,"%.15f\t%.15f\t%.15f\n",t0,y0,alpha);

    for (i=1;i<=N;i++) 
    { 
        k1=dt*f(t,y,u);
        l1=dt*u;
        k2=dt*f(t+0.5*dt,y+0.5*l1,u+0.5*k1);
        l2=dt*(u+0.5*k1);
        k3=dt*f(t+0.5*dt,y+0.5*l2,u+0.5*k2);
        l3=dt*(u+0.5*k2);
        k4=dt*f(t+dt,y+l3,u+k3);
        l4=dt*(u+k3);
        u=u+(k1+2.0*k2+2.0*k3+k4)/6.0; 
        y=y+(l1+2.0*l2+2.0*l3+l4)/6.0; 
        t=t0+dt*i;

      fprintf(output,"%.15f\t%.15f\t%.15f\n",t,y,alpha);

    } 
    fprintf(output,"\n");
    yy1=y-y1;
    t=t0;
    y=y0;
    alpha2=alpha+dalpha;
    u=alpha2;
    error=1e16;
    while(fabs(error)>1e-10){
    fprintf(output,"%.15f\t%.15f\t%.15f\n",t0,y0,alpha);

    for (i=1;i<=N;i++) 
    { 
        k1=dt*f(t,y,u);
        l1=dt*u;
        k2=dt*f(t+0.5*dt,y+0.5*l1,u+0.5*k1);
        l2=dt*(u+0.5*k1);
        k3=dt*f(t+0.5*dt,y+0.5*l2,u+0.5*k2);
        l3=dt*(u+0.5*k2);
        k4=dt*f(t+dt,y+l3,u+k3);
        l4=dt*(u+k3);
        u=u+(k1+2.0*k2+2.0*k3+k4)/6.0; 
        y=y+(l1+2.0*l2+2.0*l3+l4)/6.0; 
        t=t0+dt*i;

      fprintf(output,"%.15f\t%.15f\t%.15f\n",t,y,alpha);

    } 
    fprintf(output,"\n");
    yy2=y-y1;

    alpha = alpha2-yy2*(alpha2-alpha1)/(yy2-yy1); 
    alpha1=alpha2;
    alpha2=alpha;
    error=alpha2-alpha1;
    printf("alpha=%.15f,\terror=%.15f\n",alpha,error);
    t=t0;
    y=y0;
    u=alpha2;
    yy1=yy2;
}
    
     fclose(output);//closing file.

} 



// Define function

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

fun=-dy*dy-y+log(t);
    return fun;
} 


double analyticalsol(double x) 
{ 
    double y;

y=log(x);
    return y;
} 

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.

}