% Example 3.3
time = [ -0.2400 -0.0600  0.1200  0.3000  0.4800  0.6600  0.8400 ...
          1.0200  1.2000  1.3800  1.5600  1.7400 ]';
Y005 = [ 0.000424  -0.001633  0.004100  0.032023  0.101828  0.190799...
        0.279839  0.332546  0.356659  0.357065  0.357608  0.359628 ]';
Y01 = [-0.004206  -0.004206  0.005432  0.031767  0.102075  0.190454...
        0.277610  0.338322  0.358413  0.355197  0.359953  0.362502 ]';
Y05 = [ -0.004134  -0.003218  0.031233  -0.000487  0.089262  0.229876...
         0.282059  0.358537  0.340391  0.353100  0.342752  0.329812 ]'

dt = time(2) - time(1);
time_fine = [ -.24:0.006:1.8 ]';
qhat005 = Stolz( dt, x, Y005 );
qhat01 = Stolz( dt, x, Y01 );
qhat05 = Stolz( dt, x, Y05 );
qexact = q_triangle(time - dt/2);
figure(96); clf; 
plot( time - dt/2, qhat005, 'ro', ...
      time - dt/2, qhat01, 'bd', ...
      time - dt/2, qhat05, 'ks', ...
      time_fine, q_triangle(time_fine), ... 
      'k--', 'linewidth', 1.0)
legend('\sigma = 0.002', '\sigma = 0.004', '\sigma = 0.018', 'Exact' )
xlabel('time')
ylabel('Heat Flux')
for i = 1:length(time)
    fprintf(' %5.2f     %9.6f     %9.6f     %9.6f\n', ...
             time(i)-dt/2, qhat005(i), qhat01(i), qhat05(i) );
end

% fig=figure(96);
% print(fig,fullfile(pwd, 'Figure3_4.tif'),'-dtiff','-r600')