// to approximate the solution of nonlinear boundary value problem with finite difference
// y"=f(x,y,y') , a<x<b , y(a)=alfa	  , y(b)=beta;
// Numerical Analysis -- Richard L. Burden -- J.Douglas Faires
// Page 683 9th edition
 


class f
{
	
    public double func(double x,double y,double ydot)	
	{
	return 	(32.0+(2.0*x*x*x)-y*ydot)/8.0;
	}
}

class fy
{
    public double func(double x,double y,double ydot)
    {
	return (-ydot)/8.0;
	}	
}

class fydot
{
    public double func(double x,double y,double ydot)
    {
	return (-y)/8.0;    
	} 	
}

public class nonlinear_shoting_with_finite_difference
{
f ff;
fy ffy;
fydot ffydot;
double a;
double b;
double alfa;
double beta;
int N;
                 
	
  public nonlinear_shoting_with_finite_difference(double aa,double bb,double aalfa,double bbeta, int NN,f fff,fy fffy,fydot fffydot)
  {
  a=aa;
  b=bb;
  alfa=aalfa;  
  beta=bbeta;
  ff=fff;
  ffy=fffy;
  ffydot=fffydot;
  N=NN;
  }              
   
  public  double norm(double v[])
  {
   // vector norm
   double total=0;;
   for(int i=0;i<v.length;i++)
    {total+=v[i]*v[i];}
   return Math.sqrt(total);
  }             
  public double exact(double x)
  {
  return x*x+(16.0/x); 	  
  }
  
  
     public void solution()
     {
	 double[] w=new double [N+2];    
	 double[] aa=new double[N+1];
	 double[] bb=new double[N+1];
	 double[] dd=new double[N+1];
	 double[] cc=new double[N+1];
	 double[] xx=new double[N+2];
	 double[] l =new double[N+1];
	 double[] z =new double[N+1];
	 double[] u =new double[N+1];
	 double[] v =new double[N+1];
	 double x=0;
	 double h=(b-a)/(N+1);    
	 w[0]=alfa;
	// System.out.println(w[0]);
	 w[N+1]=beta;
	 double t=0;    
	     for(int i=1;i<=N;i++)
	     {
		 w[i]=alfa+(i*((beta-alfa)/(b-a))*h);    
		 }
	     int k=1;
		      while(k<1000)
		      {
			  x=a+h;
			  t=(w[2]-alfa)/(2.0*h);
			  aa[1]=2.0+h*h*ffy.func(x,w[1],t);
			  bb[1]=-1.0+(h/2.0)*ffydot.func(x,w[1],t);
			  dd[1]=-(2.0*w[1]-w[2]-alfa+h*h*ff.func(x,w[1],t));
			  
			      for(int i=2;i<=N-1;i++)
			      {
				  x=a+(i)*h;
				  t=(w[i+1]-w[i-1])/(2.0*h);
				  aa[i]=2.0+h*h*ffy.func(x,w[i],t);
				  bb[i]=-1.0+(h/2.0)*ffydot.func(x,w[i],t);    
				  cc[i]=-1.0-(h/2.0)*ffydot.func(x,w[i],t);      
				  dd[i]=-(2.0*w[i]-w[i+1]-w[i-1]+h*h*ff.func(x,w[i],t));    
				  }        
			      
			      x=b-h;
			      t=(beta-w[N-1])/(2.0*h);
			      aa[N]=2.0+h*h*ffy.func(x,w[N],t);
			      cc[N]=-1.0-(h/2.0)*ffydot.func(x,w[N],t);
			      dd[N]=-(2.0*w[N]-w[N-1]-beta+h*h*ff.func(x,w[N],t));
			      
			      l[1]=aa[1];
			      u[1]=bb[1]/aa[1];
			      z[1]=dd[1]/l[1];
			      
			      for(int i=2;i<=N-1;i++)
			      {
				  l[i]=aa[i]-cc[i]*u[i-1];    
				  u[i]=bb[i]/l[i];
				  z[i]=(dd[i]-cc[i]*z[i-1])/l[i];    
				  }
			      l[N]=aa[N]-cc[N]*u[N-1];
			      z[N]=(dd[N]-cc[N]*z[N-1])/l[N];
			      
			      v[N]=z[N];
			      w[N]+=v[N];
			      
			      for(int i=N-1;i>=1;i--)
			      {
				     v[i]=z[i]-u[i]*v[i+1];
			         w[i]+=v[i];
				  }
			      
			      if(norm(v)<0.1){break;} 
			      k++;
			  }
	        double[] depo=new double[N+2];
	  		for(int i=0;i<=N+1;i++)
	  		{xx[i]=a+i*h;}    
	     
	  		for(int i=0;i<=N+1;i++)
	  		{depo[i]=exact(xx[i]);
		  		System.out.println(xx[i]+ " "+w[i]+" "+exact(xx[i])+" "+Math.abs(exact(xx[i])-w[i]));}
	         Plot pp=new Plot(xx,w);
	         pp.addData(xx,depo);
	         pp.plot();
	  		
	 }             
              
                 
     public static void main(String args[])
     {
	 double aa=1.0;
	 double aalfa=17.0; 
	 double bb=3.0;
	 double bbeta=43.0/3.0;
	 int NN=200;
	 f fff=new f();
	 fy fffy=new fy() ;
	 fydot fffydot=new fydot() ;    
	 
	     
	 nonlinear_shoting_with_finite_difference fd=new nonlinear_shoting_with_finite_difference(aa,bb,aalfa,bbeta,NN,fff,fffy,fffydot);   
	 fd.solution();   
	     
	     
	 }         
                  
                     
                    
  







}