/* Solution of y''=y, y[0]=0.0, y[1]=1.0 by shooting
 Since this is a linear problem we can use two independent 
 IVP solutions to determine the solution 

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

The exact solution is y(x)=sin(x)/sin(1.)
*/

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

#define pi acos(-1.0L)
#define sqr(x) ((x)*(x))
#define m 1

#define h 0.001
#define M 1000
#define FILE1 "sol1.txt"
#define FILE2 "sol2.txt"

double F(double a, double b){
  return(-sqr(m*pi)*b);
}

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

void RK4(double *X, double *Y){
/*    double a, b;

    a = F(*X,*Y);
    b = G(*X,*Y);

    (*X) = (*X) + a;
    (*Y) = (*Y) + b;*/

    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]);

}

main(){
    double X, Y, X1, Y1, X2, Y2; 
    double Yi, Yi1, Yi2, Yf, Yf1, Yf2; 
    double c1, c2;

    double sol_exact, error;

    int i;
    FILE *out_file1;
    out_file1=fopen(FILE1, "w");
          
// The boundary value problem: Yi=0.0, Yf=1.0

// I store the value of Y at the boundary t=0

    Yi=0.0;

// I store the value of Y at the boundary t=1

    Yf=1.0; 

// First Initial Value Problem 
    X=1.0; Y=1.0;
// I give the boundary condition for Y at t=0 and provide a value for the slope at that point (here X=1.0)
    
// I store the value of Y at the boundary t=0 for the first problem

    Yi1=Y;

    for(i=1;i<=M;i++){
	RK4(&X, &Y);
	fprintf(out_file1, "%18.16g %18.16g\n",X, Y);
    }

// I store the (computed) value of Y at the boundary t=1 for the first problem

    Yf1 = Y;

    fclose(out_file1);	

    FILE *out_file2; 
    out_file2=fopen(FILE2, "w");


// Second Initial Value Problem
    X=-1.0; Y=1.0; 

// I store the value of Y at the boundary t=0 for the second problem

    Yi2 = Y; 

// I give the same boundary condition for Y at t=0 and provide another value for the slope at that point (here X=-1.0) I thus obtain a linearly indep. solution

    for(i=1;i<=M;i++){
	RK4(&X, &Y);
	fprintf(out_file2, "%18.16g %18.16g\n",X, Y);
    }

// I store the (computed) value of Y at the boundary t=1 for the second problem

    Yf2 = Y; 

    fclose(out_file2);

/* I know have two solutions to the problem X'(t)=F(X(t),Y(t)) and 
   Y'(t)=G(X(t),Y(t))   on the interval 0<= t <= 1 . 
I use the two solutions at t=1 to compute the constants c1 and c2 and express
the solution as c1*y1(t)+c2*y2(t) */

    c1 = 1.0;
    c2 = -Yf1/Yf2;

// at t=0
    i=0;
    sol_exact = sin((double)(m*pi*i*h));
    printf("%g %g %g\n", i*h, c1*Yi1+c2*Yi2,fabs(sol_exact-c1*Yi1-c2*Yi2)); 
// at t>0

    out_file1=fopen(FILE1, "r");
    out_file2=fopen(FILE2, "r");

    for(i=1;i<=M;i++){
      sol_exact = (2.0/pi)*sin((double)(m*pi*i*h));

	fscanf(out_file1, "%lf %lf", &X1, &Y1);
	fscanf(out_file2, "%lf %lf", &X2, &Y2);
	
		error = fabs(sol_exact - (c1*Y1 + c2*Y2) );
	
		printf("%g\t%g\t%18.16g\n", i*h, c1*Y1+c2*Y2,error);

		//	printf("%g\t%g\n",i*h, c1*Y1+c2*Y2);
    }

    fclose(out_file1);
    fclose(out_file1);



}