//-------------------
// Solve a linear BVP
//
// -x''(t)-(1+t^2)x(t)=1,    -1<t<1
//
// x(-1)=x(1)=0
//
// with the finite difference method.
//--------------------
/* 
compile with	
gcc -o FinDiff FinDiff.c  -L/usr/local/dislin -ldislnc -lm -lX11 -lgsl -lgslcblas
*/

#include "/usr/local/dislin/dislin.h"
#include <stdio.h>
#include <math.h>
#include </usr/include/gsl/gsl_matrix.h>
#include </usr/include/gsl/gsl_vector.h>
#include </usr/include/gsl/gsl_linalg.h>

//----------------------------------------------------------
int i;
// Timeinterval [-L, L], L=1;
double L = 1.0;
// Amount of points on [-L, L]
double M = 200.0;

double h=0.01;

double fun_q(int i){
  double t = -L + i * h;
  return -(1.0+pow(t,2));
}
//-------------------------------------------------------------
main ()
{
// Calculate the step size
double h=(2.0*L)/M;
// define the system A*y=b, where A is tridiagonal

gsl_vector * x = gsl_vector_alloc (M+1); // allocate vector x
gsl_vector * y = gsl_vector_alloc (M-1); // allocate vector x
gsl_vector * b = gsl_vector_alloc (M-1); // allocate vector b
gsl_vector * d = gsl_vector_alloc (M-1); // main diagonal of the matrix A
gsl_vector * d1 = gsl_vector_alloc (M-2); // off-diagonal vector of the matrix A

gsl_vector_set (x, 0, 2.0); // boundary condition on the left
gsl_vector_set (x, M, 1.0); // boundary condition on the right

/*for(i=0; i<M+1; i++){
printf (" X(%d) = %g\n", i, 
                   gsl_vector_get(x, i));
}
*/
// define vector b with the elements b_i=-h^2, i=0,M
gsl_vector_set_all (b, -pow(h,2));

// define the main diagonal of the matrix A
for (i=0; i<M-1; i++){
 gsl_vector_set (d, i, -(2.0+pow(h,2)*fun_q(i)));
}

// define the off-diagonal of symmetric A 
gsl_vector_set_all (d1, 1.0);

gsl_linalg_solve_symm_tridiag (d, d1, b,y);

for (i=1; i<M; i++){
gsl_vector_set (x,i,gsl_vector_get(y, i-1));
}
// Plot the result

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

  disini ();
  pagera ();
  hwfont ();
setrgb(0,0,0);
axspos(450,1800);
axslen(1800,1200);
name("t","x");
name("x(t)","y");
graf(-L-0.4,L+0.4,-L-0.4,0.2,0.0,1.2,0.0,0.2);
color("blue");
hsymbl(6);
      for (i=0; i<M; i++) {
        rlsymb(21, -L + i * h, gsl_vector_get(x, i));
      }
//endgrf();

 gsl_vector_free (x);gsl_vector_free (y);gsl_vector_free (b);gsl_vector_free (d);gsl_vector_free (d1);
disfin ();
}