% Example 3.4
time = [ -.24: 0.06: 1.8]';
x = 0.5;
  
Y = [ ...
0.0000000	0.0000000 0.0000000 0.0000000 0.0000000 0.0003952 0.0034040 ...
0.0099142 0.0200152 0.0337153 0.0510153 0.0719153 0.0964153 0.1245153 ...
0.1562153 0.1907249 0.2236072 0.2530868 0.2789849 0.3012847 0.3199847 ...
0.3350847 0.3465847 0.3544847 0.3587847 0.3598799 0.3599888 0.3599990 ...
0.3599999 0.3600000 0.3600000 0.3600000 0.3600000 0.3600000 0.3600000 ...
]';


dt = time(2) - time(1);
time_fine = [ -.24:0.006:1.8 ]';
qhat = Stolz( dt, x, Y );
qexact = q_triangle(time - dt/2);
figure(96); clf; 
plot( time - dt/2, qhat, 'ro-', ...
      time_fine, q_triangle(time_fine), ... 
      'k--', 'linewidth', 1.0)

xlabel('time')
ylabel('Heat Flux')
axis( [ -.5 2 -.4 .9 ] )
for i = 1:length(time)
    fprintf(' %5.2f     %9.6f\n', ...
             time(i)-dt/2, qhat(i) );
end

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