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

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,r=k/h/h;
	double La[M],Lb[M-1],Lc[M-1],F[M],Y[M],Z[M];
	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;
		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;
		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);
		tdmult(M-1,&La[1],&Lb[1],&Lc[1],&Y[1],&Z[1]);
		for(int m=1;m<M;m++){
			Y[m]+=r*(Z[m]+F[m]);
		}
	}
	return 0;
}