// Linear Shooting Method

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

void ode2rk4(double (*ddfunc)(double x, double y, double yd),
             double x0, double y0, double yp0, double xe, int n, double y[]);
double ddfunc1(double x, double y, double yd);
double ddfunc2(double x, double y, double yd);
double realfunc(double x);

void main()
{
   int i, n=20;  // n은 분할 수
   double a=1, alpha=-1.0/2.0, b=2, beta=log(2.0);    // 초기치
   double y1[100], y2[100], result[100];
   double x, h=(b-a)/n;

   ode2rk4(ddfunc1, a, alpha, 0, b, n, y1);
   ode2rk4(ddfunc2, a, 0, 1, b, n, y2);

   for(i=0; i<=n; i++)
      result[i] = y1[i] + (beta-y1[n])/y2[n] * y2[i];
   
   printf(" Xi      y1(1,i)      y2(1,i)        Wi          y(Xi)      |y(Xi)-Wi|\n");   

   for(i=0; i<=n; i++) {
      x = a + h*i;
      printf("%4.2f  %11.9f  %11.9f  %11.9f  %11.9f  %11.9e\n",
         x, y1[i], y2[i], result[i], realfunc(x), fabs(realfunc(x)-result[i]) );
   }
}

// 정확한 함수
double realfunc(double x)
{
   return (exp(2)/(exp(4)-1.0) * (exp(2*x)-exp(-2*x))+x);
}

// 입력 함수 부분 1
// y'' = p(x)y' + q(x)y + r(x)
double ddfunc1(double x, double y, double yd)
{
   return ( -4.0/x * yd - 2.0/(x*x) * y -2.0*log(x)/(x*x) );
}

// 입력 함수 부분 2
// y'' = p(x)y' + q(x)y
double ddfunc2(double x, double y, double yd)
{
   return ( -4.0/x * yd - 2.0/(x*x) * y );
}

// 2계 상미분 방정식의 해를 구하는 함수 - 4차 RK 법 이용
// 초기치 : y(x0)=y0, y'(x0)=yp0
// 구간의 끝점 : xe
// 분할수 : n
// 결과는 y[]에 저장 : y[0]=y0, y[1]=y1, ....    
void ode2rk4(double (*ddfunc)(double x, double y, double yd),
             double x0, double y0, double yp0, double xe, int n, double y[])
{
   int i;
   double  h, k, xi, ui, udi, a, b, c, d, e, f, beta, delta;

   h = (xe-x0)/n;
   k = h*0.5;
   y[0] = ui = y0;
   udi = yp0;
   for(i=1, xi=x0; i<=n; i++, xi+=h) {
      a = (*ddfunc)(xi, ui, udi)*k;
      beta = (udi+a/2)*k;
      b = (*ddfunc)(xi+k, ui+beta, udi+a)*k;
      c = (*ddfunc)(xi+k, ui+beta, udi+b)*k;
      delta = h*(udi+c);
      d = (*ddfunc)(xi+h, ui+delta, udi+2*c)*k;
      e = (a+b+c)/3.;
      f = (b+c+d)/3.;
      ui += (udi+e)*h;
      udi += (e+f);
      y[i] = ui;
   }
}