import mt.BandMatrix;
import mt.DenseVector;

/**
 * @author sg
 */
public class LinearFDM  {
	
	public double[] solve(Function1 p, Function1 q, Function1 r, double xa, double ya, double xb, double yb, int n){
		double[] result=new double[n+1];
		result[0]=ya;
		result[n]=yb;

		double h=(xb-xa)/n;
		double hh=h/2;
		double h2=h*h;
		
		double[] x=new double[n+1];
		x[0]=xa;
		x[n]=xb;
		for(int i=1; i<n; i++){
			x[i]=x[i-1]+h;
		}
		
		BandMatrix m=new BandMatrix(n-1, 1, 1);
		m.set(0, 0, 2+h2*q.f(x[1]));
		for(int i=1; i<n-1; i++){
			m.set(i, i, 2+h2*q.f(x[i+1]));
			m.set(i, i-1, -1-hh*p.f(x[i+1]));
			m.set(i-1, i, -1+hh*p.f(x[i]));	
		}
		
		DenseVector b=new DenseVector(n-1);
		b.set(0, -h2*r.f(x[1])+(1+hh*p.f(x[1]))*ya );
		for(int i=1; i<n-2; i++){
			b.set(i, -h2*r.f(x[i+1]));
		}
		b.set(n-2, -h2*r.f(x[n-1])+(1-hh*p.f(x[n-1]))*yb );
		
		DenseVector sol=new DenseVector(n-1);
		m.solve(b, sol);
		System.arraycopy(sol.getData(), 0, result, 1,  n-1);
		return result;
	}
}