function printTimeEvolution(timeVector, dataMatrix, ...
    dataMatrixReference, labelList)
% compare two solutions after approximation
%
% input data as constructed by readResultData

%% refine data (if multiple points per time)
uniqueTimeVector = unique(timeVector);
if length(timeVector) > length(uniqueTimeVector)
    newDataMatrix = zeros(size(dataMatrix, 1), length(uniqueTimeVector));
    newDataMatrixReference = newDataMatrix;
    for iTime = 1:length(uniqueTimeVector)
        thisFilter = timeVector == uniqueTimeVector(iTime);
        newDataMatrix(:, iTime) = mean(dataMatrix(:, thisFilter), 2);
        newDataMatrixReference(:, iTime) = ...
            mean(dataMatrixReference(:, thisFilter), 2);
    end
    timeVector = uniqueTimeVector;
    dataMatrix = newDataMatrix;
    dataMatrixReference = newDataMatrixReference;
end

figure
subplot(2,1,1)
hold on

colorMatrix = [0 0 0; 1 0 0; 0 1 0; 0 0 1; ...
    1 1 0; 1 0 1; 0 1 1; 0.5 1 0; 1 0.5 0; 0 0.5 1; 0 1 0.5; 0.5 0 1; 1 0 0.5];

% reference
for iData = 1:size(dataMatrix, 1)
    plot(timeVector, dataMatrixReference(iData, :), 'r-', ...
        'color', colorMatrix(iData,:), 'lineWidth', 3);
end
if exist('labelList', 'var')
    legend(labelList)
end

%  simulation
for iData = 1:size(dataMatrix, 1)
    plot(timeVector, dataMatrix(iData, :), 'r-', ...
        'color', colorMatrix(iData,:), 'lineWidth', 1)
end
xlim([min(timeVector) max(timeVector)])

% error
subplot(2,1,2)
hold on

for iData = 1:size(dataMatrix, 1)
    thisError = dataMatrix(iData, :)./dataMatrixReference(iData, :) - 1;
    plot(timeVector, thisError, 'r--', ...
        'color', colorMatrix(iData,:), 'lineWidth', 1);
end
if exist('labelList', 'var')
    legend(labelList)
end
xlim([min(timeVector) max(timeVector)])