function la9()
	yinit  = [  1;  1; 1; 1; 1 ];
	lambda = [ -5; -1; 0; 1; 2 ];
	f = @(t, y) lambda .* y;
	h = 0.1;
	tspan = [ 0, 2 ];
	
	exact = exp(lambda * tspan(2)) .* yinit;
	
	[ T, Y ] = ExplicitRungeKutta(f, tspan, yinit, h);
	
	error_abs_rk4 = Y(:, end) - exact;
	error_rel_rk4 = (Y(:, end) - exact) ./ abs(exact);
	
	subplot(1, 2, 1);
	plot(T, Y);
	axis square;
	
	[ T, Y ] = Explicit3o8R(f, tspan, yinit, h);
	
	error_abs_3o8 = Y(:, end) - exact;
	error_rel_3o8 = (Y(:, end) - exact) ./ abs(exact);
	
	subplot(1, 2, 2);
	plot(T, Y);
	axis square;
	
	print la9.eps;
	
	stab = @(x) [ '  ' x '    ' ];
	error = [ error_abs_rk4, error_rel_rk4, error_abs_3o8, error_rel_3o8 ];
	errNames = { 'RK4-abs', 'RK4-rel', '3o8-abs', '3o8-rel' };
	
	fprintf(' Verf/lambda %s\n', tab(lambda));
	fprintf('--------------+---------------+---------------+---------------+---------------+---------------\n');
	for k = 1:length(errNames)
		fprintf('%s%s\n', stab(errNames{k}), tab(error(:, k)));
	end
end

function [ T, Y ] = ExplicitRungeKutta(f, tspan, yinit, h)
	N = round(diff(tspan) / h);
	T = linspace(tspan(1), tspan(2), N+1);
	h = diff(tspan) / N;
	Y = zeros(length(yinit), length(T));
	Y(:, 1) = yinit;
	
	t    = tspan(1);
	y    = yinit;
	Tend = tspan(2);
	
	c1 = 0;
	c2 = 1/2;
	c3 = 1/2;
	c4 = 1;
	
	a21 = 1/2;
	a32 = 1/2;
	a43 = 1;
	
	b = [ 1/6 1/3 1/3 1/6 ];
	
	for k = 2:length(T)
		t = T(k);
		
		Y1 = y;
		k1 = f(t + c1 * h, Y1);
		
		Y2 = y + h * a21 * k1;
		k2 = f(t + c2 * h, Y2);
		
		Y3 = y + h * a32 * k2;
		k3 = f(t + c3 * h, Y3);
		
		Y4 = y + h * a43 * k3;
		k4 = f(t + c4 * h, Y4);
		
		K = [ k1, k2, k3, k4 ];
		
		y = y + h * K * b';
		
		Y(:, k) = y;
	end
end

function [ T, Y ] = Explicit3o8R(f, tspan, yinit, h)
	N = round(diff(tspan) / h);
	T = linspace(tspan(1), tspan(2), N+1);
	h = diff(tspan) / N;
	Y = zeros(length(yinit), length(T));
	Y(:, 1) = yinit;
	
	t    = tspan(1);
	y    = yinit;
	Tend = tspan(2);
	
	c1 = 0;
	c2 = 1/3;
	c3 = 2/3;
	c4 = 1;
	
	a21 =  1/3;
	a31 = -1/3; a32 =  1;
	a41 =  1;   a42 = -1; a43 = 1;
	
	b = [ 1/8 3/8 3/8 1/8 ];
	
	for k = 2:length(T)
		t = T(k);
		
		Y1 = y;
		k1 = f(t + c1 * h, Y1);
		
		Y2 = y + h * a21 * k1;
		k2 = f(t + c2 * h, Y2);
		
		Y3 = y + h * (a31 * k1 + a32 * k2);
		k3 = f(t + c3 * h, Y3);
		
		Y4 = y + h * (a41 * k1 + a42 * k2 + a43 * k3);
		k4 = f(t + c4 * h, Y4);
		
		K = [ k1, k2, k3, k4 ];
		
		y = y + h * K * b';
		
		Y(:, k) = y;
	end
end

function res = tab(vec)
	res = '';
	for k = 1:length(vec)
		res = [ res ' |' sprintf('%14g', vec(k)) ];
	end
end
