/* Solution of y''=-f(y), y[0]=0.0, y[1]=1.0 by shooting
v1.0 For a linear problem we can use two independent 
 IVP solutions to determine the solution 

v2.0 To test non-linear shooting, I will now solve this problem 
     by imposing the boundary condition at t=0 and solving for the slope, i.e. 
     the value of x(0).

v2.1 I replace f(y)=y by a non-linear function of y

 I write the problem as y'=x, x'=-f(y) with F(X,Y)=-f(y), G(X,Y)=x 
 So x(t) is the slope of y(t) */

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

#define h 0.00001
#define M 100000
#define TOL 1.0e-07

#include "f.h" // This file contains the specification for the function f(y)
    
double F(double a, double b){
    return(-f(b));
}

double G(double a, double b){
    return(a);
}

void RK4(double *X, double *Y){
    double temp_X, temp_Y, temp_K[4], temp_L[4];

    temp_K[0]=h*F(*X, *Y);
    temp_L[0]=h*G(*X, *Y);

    temp_X=(*X)+0.5*temp_K[0];
    temp_Y=(*Y)+0.5*temp_L[0];

    temp_K[1]=h*F(temp_X, temp_Y);
    temp_L[1]=h*G(temp_X, temp_Y);


    temp_X = (*X) + 0.5*temp_K[1];
    temp_Y=  (*Y) + 0.5*temp_L[1];

    temp_K[2]=h*F(temp_X, temp_Y);
    temp_L[2]=h*G(temp_X, temp_Y);

    temp_X = (*X) + temp_K[2];
    temp_Y=  (*Y) + temp_L[2];

    temp_K[3]=h*F(temp_X, temp_Y);
    temp_L[3]=h*G(temp_X, temp_Y);


    (*X) = (*X) + (1./6.)*(temp_K[0]+2.0*(temp_K[1]+temp_K[2]) + temp_K[3]);
    (*Y) = (*Y) + (1./6.)*(temp_L[0]+2.0*(temp_L[1]+temp_L[2]) + temp_L[3]);

}

// This function finds the value of the solution (X, Y) at the endpoint t=1

void sol(double *X, double *Y){
    int i;
    for(i=1;i<=M;i++){
	RK4(X, Y);
    }
}





main(){
    double X, Y; 
    double Xi, Xi1, Xin, dx, Yi, Yf, Yf1, Yf2, Yfn;
    int k;

// impose the boundary condition on Y at t=0
    Yi=0.0; 
    Y=Yi; 
    
// the boundary condition at t=1
    Yf=1.0; 


// This part of the code computes the slope at t=0==========================

// choose initial slope and increment
    Xi=1.0; dx=0.1; 
//compute corresponding boundary value at t=1
    X=Xi; Y=Yi; sol(&X,&Y); Yf1=Y;

//compute boundary value at t=1 for slope Xi+dx
    Xi1=Xi+dx; X=Xi1; Y=Yi; sol(&X, &Y); Yf2=Y;

//Now I can set up the Newton-Raphson iteration
    while( (fabs(Yf2-Yf)>TOL) || ((Xi1-Xi)>TOL) ){
    Xin = Xi1 - ((Xi1-Xi)*(Yf2-Yf)/(Yf2-Yf1));
    X=Xin; Y=Yi; sol(&X, &Y); 
    Yf1=Yf2; Yf2=Y; Xi=Xi1; Xi1=Xin; 
//    printf("%g %g\n",Xin, Y);
}

//============================================================================
    
// Now I know the slope at t=0 so I can compute the solution from t=0 to t=1
    X=Xin; Y=Yi; 
    for(k=1;k<=M;k++){
	RK4(&X,&Y);
	printf("%g %g %g\n",k*h, X, Y);
    }

    

}