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

void tridiaglu(int n,double a[n],
	double b[n-1],double c[n-1]){
	for(int k=1;k<n;k++){
		b[k-1]=b[k-1]/a[k-1];
		a[k]=a[k]-b[k-1]*c[k-1];
	}
}
void lusolve(int n,double a[n],
	double b[n-1],double c[n-1],
	double f[n],double y[n]){
	double z[n];
	z[0]=f[0];
	for(int k=1;k<n;k++) z[k]=f[k]-b[k-1]*z[k-1];
	y[n-1]=z[n-1]/a[n-1];
	for(int k=n-2;k>=0;k--) y[k]=(z[k]-c[k]*y[k+1])/a[k];
}
void tdmult(int n,double a[n],
	double b[n-1],double c[n-1],double x[n],double y[n]){
	y[0]=a[0]*x[0]+c[0]*x[1];
	for(int k=1;k<n-1;k++){
		y[k]=b[k-1]*x[k-1]+a[k]*x[k]+c[k]*x[k+1];
	}
	y[n-1]=b[n-2]*x[n-2]+a[n-1]*x[n-1];
}
double v(double t,double x){
	return 0.01;
}
double p(double t,double x){
	return 0;
}
double q(double t,double x){
	return 0;
}
double f(double t,double x){
	return 0;
}
double yo(double x){
	if(x<1) return 2*x;
	else return 3-x;
}
double A(double t){
	return 0;
}
double B(double t){
	return 1;
}

int main(){
	int M=100,N=10000;
	double a=0,b=2,h=(b-a)/M;
	double t0=0,T=4,k=(T-t0)/N,r2=k/h/h/2;
	double La[M],Lb[M-1],Lc[M-1],F[M],Y[M],Z[M];
	double ImLa[M],ImLb[M-1],ImLc[M-1];
	for(int m=1;m<M;m++){
		double xm=a+m*h;
		Y[m]=yo(xm);
	}
	for(int n=0;;n++){
		double tn=t0+n*k;
		for(int m=1;m<M;m++){
			double xm=a+m*h;
			La[m]=-2*v(tn,xm)+h*h*q(tn,xm);
			Lb[m-1]=v(tn,xm)-h*p(tn,xm)/2; 
			Lc[m]=v(tn,xm)+h*p(tn,xm)/2;
		}
		F[1]=f(tn,a)+Lb[0]*A(tn);
		for(int m=1;m<M-1;m++){
			double xm=a+m*h;
			F[m]=f(tn,xm);
		}
		F[M-1]=h*h*f(tn,b)+Lc[M-1]*B(tn);
		if(n>0){
			for(int m=1;m<M;m++){
				ImLa[m]=1-r2*La[m];
				ImLb[m-1]=-r2*Lb[m-1];
				ImLc[m]=-r2*Lc[m];
				Z[m]+=r2*F[m];
			}
			tridiaglu(M-1,&ImLa[1],&ImLb[1],&ImLc[1]);
			lusolve(M-1,&ImLa[1],&ImLb[1],&ImLc[1],&Z[1],&Y[1]);
		}
		if(n%100==0){
			printf("# tn=%g\n",tn);
			printf("%g %g\n",a,A(tn));
			for(int m=1;m<M;m++){
				double xm=a+m*h;
				printf("%g %g\n",xm,Y[m]);
			}
			printf("%g %g\n\n\n",b,B(tn));
		}
		if(n>=N) break;
		tdmult(M-1,&La[1],&Lb[1],&Lc[1],&Y[1],&Z[1]);
		for(int m=1;m<M;m++){
			Z[m]=Y[m]+r2*(Z[m]+F[m]);
		}
	}
	return 0;
}