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

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

#define max 2000


int max_steps = 4.;   // time interval [0, max_steps]
int i, j, nn;                    /* Loop counter                   */
    float T[max];    
    float X1[max];              
    float Y1[max];
    float X2[max];               
    float Y2[max];
    float Xf[max];
double h = 0.005;
// initial conditions
double x = 1.25; // boundary condition on the left side of the interval [0, 4]
double y = 0.0; 
double x2 = 0.0;
double y2 = 1.0;
double b= -0.95; // boundary condition on the right side of the interval [0, 4]


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

// r.h.s. for the function u
double funct_x1(double t, double x, double y) {
  return y;
}

double funct_y1(double t, double x, double y) {
  return (2.0*t/(1.+pow(t,2)))*y-(2.0/(1.+pow(t,2)))*x+1.0;
}

// r.h.s. for the function v
double funct_x2(double t, double x2, double y2) {
  return y2;
}

double funct_y2(double t, double x2, double y2) {
  return (2.0*t/(1.+pow(t,2)))*y2-(2.0/(1.+pow(t,2)))*x2;
}
// RK4 method for the equation for u(t)
void rk4u(double h, double max_steps, double t, double x, double y)
{
    T[0] = t;
    X1[0] = x;
    Y1[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  = X1[j];
      y  = Y1[j];

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

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

    k4x = h * funct_x1(t+h, x + k3x, y + k3y);
    k4y = h * funct_y1(t+h, x + k3x, y + k3y);

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

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

    T[j+1] = t + h;
    j++;
   }
}

// RK4 method for the equation for v(t)
void rk4v(double h, double max_steps, double t, double x2, double y2)
{
    T[0] = t;
    X2[0] = x2;
    Y2[0] = y2;
    j    = 0;
while(T[j] < max_steps)
    {
      if( (T[j] + h) > max_steps ) h = max_steps - T[j];   
      t  = T[j];
      x2  = X2[j];
      y2  = Y2[j];

    k1x = h * funct_x2(t, x2, y2);
    k1y = h * funct_y2(t, x2, y2);
    
    k2x = h * funct_x2(t+0.5*h, x2 + 0.5 * k1x, y2 + 0.5 * k1y);
    k2y = h * funct_y2(t+0.5*h, x2 + 0.5 * k1x, y2 + 0.5 * k1y);
    

    k3x = h * funct_x2(t+0.5*h, x2 + 0.5 * k2x, y2 + 0.5 * k2y);
    k3y = h * funct_y2(t+0.5*h, x2 + 0.5 * k2x, y2 + 0.5 * k2y);

    k4x = h * funct_x2(t+h, x2 + k3x, y2 + k3y);
    k4y = h * funct_y2(t+h, x2 + k3x, y2 + k3y);

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

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

    T[j+1] = t + h;
    j++;
   }
}

//---------------------------------------------------------------------

main ()
{ 
// calculate u(t)
rk4u(h,max_steps,t,x,y);
t=0.0;
//calculate v(t)
rk4v(h,max_steps,t,x2,y2);
nn=j;
// claculate x(t)=u(t)+(b-u(nn))*v(t)/v(nn)
for (i = 0; i <=nn; i++) Xf[i] = X1[i] + (b - X1[nn]) * X2[i] / X2[nn];

//-----------------------------------------------------------------

// plot with DISLIN
 winsiz (600, 600);
  page (2600, 2600);
  sclmod ("full");
  scrmod ("revers");
  metafl ("xwin");

  disini();
  pagera();
  texmod ("on");
axspos(450,1800);
  axslen(1800,1200);

  name("t","x");
  name("u(t), v(t), x(t)","y");

 
  graf(0.,max_steps,0.,1.0,-max_steps,max_steps, -max_steps,1.0);
  height(50);
  title();
// plot u(t)
  color("red");
  curve(T,X1,nn);
// plot v(t)
  color("blue");
  curve(T,X2,nn);
// plot x(t)
  color("green");
  curve(T,Xf,nn);
  endgrf();
  disfin ();
}