//-------------------
// Solution of Boundary Value Problem with the shooting method.
// Equation to solve:
// x''=a*x^2 with a=1,5 on
// the interval [0,1] 
// Boundary conditions are: x(0)=4; x(1)=1;
//--------------------
/* 
compile with	
gcc -o NLinearShooting NLinearShooting.c  -L/usr/local/dislin -ldislnc -lm -lX11
*/

#include "/usr/local/dislin/dislin.h"
#include <stdio.h>
#include <math.h>

#define max 500


int max_steps = 1.;   // t in [0, max_steps]
int i, j, k, nn;                    /* Loop counter                   */

    float T[max];    
    float X[max];              
    float Y[max];
    float F[max];
    float St[max];
    float s0[max];
double h = 0.005; // time step
// initial conditions for the corresponding i.v.p.
double x = 4.0; // boundary condition on the left side of the interval [0, max_steps]
double s = -100.0; // slope to be varied
double b= 1.0; // boundary condition on the right side of the interval [0, max_steps]


double k1x, k2x, k3x, k4x, k1y, k2y, k3y, k4y;  /* Function values for RK4    */
double t=0.0;

// r.h.s. for the o.d.e's of the 1st order
double funct_x(double t, double x, double y) {
  return y;
}

double funct_y(double t, double x, double y) {
  return 1.5*pow(x,2);
}
// RK4 method
void rk4(double h, double max_steps, double t, double x, double y)
{
    T[0] = t;
    X[0] = x;
    Y[0] = y;
    j    = 0;
while(T[j] < max_steps)
    {
      if( (T[j] + h) > max_steps ) h = max_steps - T[j];   /* Calculation of the last step */
      t  = T[j];
      x  = X[j];
      y  = Y[j];

    k1x = h * funct_x(t, x, y);
    k1y = h * funct_y(t, x, y);
    
    k2x = h * funct_x(t+0.5*h, x + 0.5 * k1x, y + 0.5 * k1y);
    k2y = h * funct_y(t+0.5*h, x + 0.5 * k1x, y + 0.5 * k1y);
    

    k3x = h * funct_x(t+0.5*h, x + 0.5 * k2x, y + 0.5 * k2y);
    k3y = h * funct_y(t+0.5*h, x + 0.5 * k2x, y + 0.5 * k2y);

    k4x = h * funct_x(t+h, x + k3x, y + k3y);
    k4y = h * funct_y(t+h, x + k3x, y + k3y);

    X[j+1]= x+(k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;

    Y[j+1]= y+(k1y + 2.0 * k2y + 2.0 * k3y + k4y) / 6.0;

    T[j+1] = t + h;
    j++;
   }
}
//----------------------------------------------------------
main ()
{


  winsiz (800, 800);
  page (2600, 2600);
  sclmod ("full");
  scrmod ("revers");
  metafl ("xwin");

  disini ();
  pagera ();
  hwfont ();
  texmod ("on");

  axspos(250,2000);
  axslen(900,900);
  graf(-100.0,0.,-100.0,20.0,-10.0,80.0, -10.0,10.0);
  endgrf();

  
  axspos(1600,2000);
  axslen(900,900);
  graf(0.0,max_steps,0.0,0.2,-12.0,4.0,-12.0,2.0);
  endgrf();
// --------------------------find F(s)=x(Tend,s)-b;----------------------------
St[0]=s;
for (i=0; i<=100; i++){
t=0.0;
x=4.0;
rk4(h,max_steps,t,x,s);
nn=j;
F[i]=X[nn]-b;
St[i+1]=s+1;
s++;
}
  
// ---------------plot F(s)------------------------
  axspos(250,2000);
  axslen(900,900);
  name("s","x");
  name("F(s)","y");
  graf(-100.0,0.,-100.0,20.0,-10.0,80.0, -10.0,10.0);
  color("red");
  curve(St,F,i);
  endgrf();

// --------find zeros of the function F(s) with the Newton's method-----------------

// find initial estimation
s=-100.0;
k=0;
 for (i=0; i<=100; i++){
     if (((F[i]<0) && (F[i+1]>0)) || ((F[i]>0) && (F[i+1]<0)))
     {    s0[k]=s;
          printf("Initial guess for the Newton method: s0 =  %lf\n", s0[k]);
           printf("k = %d, s0[k] = %lf\n", k, s0[k]);
          k++;
     }
      s++;
}
//------------------------------ Newton method---------------------------------

float F1[max];
float F2[max];
float So[max];
float Sopt[max];
double ds=0.01;
k=0;
for (k=0; k<2;){
So[0]=s0[k];
for (i=0; i<=100; i++){
rk4(h,max_steps,t,x,So[i]);
nn=j;
F1[i]=X[nn]-b;
rk4(h,max_steps,t,x,So[i]+ds);
nn=j;
F2[i]=X[nn]-b;
So[i+1]=So[i]-(F1[i]*ds)/(F2[i]-F1[i]);
}
Sopt[k]=So[i];
printf("k=%d, Sopt =  %lf\n",k, Sopt[k]);
rk4(h,max_steps,t,x,Sopt[k]);
 setrgb(0,0,0);
axspos(1600,2000);
axslen(900,900);
name("t","x");
  name("x(t)","y");
graf(0.0,max_steps,0.0,0.2,-12.0,4.0,-12.0,2.0);
color("blue");
curve(T,X,nn); 
endgrf();
k++;
}
// plot solutions
   
  

  disfin ();
//----------------------------------------------------------
}