/*
author:
	seth from http://www.bierdatenbank.de
tab stops:
	2 (otherwise hardly readable!)
*/

#include "matrix.h"

void explEuler(vektor& y, matrix& S, double tau, double int_start, double int_end){
	double it=(int_end-int_start)/tau;
	for(double i=0;i<it;++i) y+=tau*(S*y).m2v();
	cout <<"tau = "<<tau<<", y = "<<y<<endl;
}

void aufgabe14(){
	cout <<"14.aufg:"<<endl;
	cout <<"========"<<endl;
	int n=5, n2=n*2;
	double k=2, m=10, schrittweite=1.e-3, int_start=0, int_end=5;
	matrix S(n2,n2);
	vektor y0(n2);
	S.zero();
	double help=k*m*(-2);
	for(int i=0;i<n;++i){
		S.set(i,i+n,1);
		S.set(i+n,i,help);
		y0.set(i,0);
		y0.set(i+n,1);
	}
	help=k*m;
	for(int i=1;i<n;++i){
		S.set(i+n-1,i,help);
		S.set(i+n,i-1,help);
	}
	for(double tau=schrittweite;tau>1.e-5;tau/=2){
		vektor y(y0);
		explEuler(y, S, tau, int_start, int_end);
		y.deallocate();
	}
	S.deallocate();
	y0.deallocate();
}

int main(){
	aufgabe14();
	return 0;
}