%
% clear previous plot
clg
%
% load output from Fortran program for 3-D plot of u(r,z) vs r and z
load h5a.out;
%
% number of r and z values
nr=21;
nz=101;
%
% r grid
r=(0:0.05:1);
%
% z grid
z=h5a(1:nz,1);
%
% read output of Fortran program
%
% dependent variable
%
% Pe = 60
for i=1:nz;
for j=1:nr;
    u60(i,j)=h5a(i,j+1);
end;
end;
%
% Pe = 120
for i=1:nz;
for j=1:nr;
    u120(i,j)=h5a(i+nz,j+1);
end;
end;
%
% Pe = 180
for i=1:nz;
for j=1:nr;
    u180(i,j)=h5a(i+2*nz,j+1);
end;
end;
%
% plot u(r,z) vs r and z
figure(1);
%
% Pe = 60
mesh(r,z,u60);
view(-135,45);
hold on;
%
% Pe = 120
mesh(r,z,u120);
view(-135,45);
hold on;
%
% Pe = 180
mesh(r,z,u180);
view(-135,45);
hold on;
%
% label three axes
xlabel(' r');
ylabel(' z');
zlabel(' u(r,z)');
title(' Graetz problem with constant wall temperature, Pe = 60, 120, 180')
%
% print output
print h5p9a.eps -deps
%
% store series solution
load h5b.out;
%
% Pe = 60
for i=1:nz;
   ua60(i,1)=h5b(i,2);
end;
%
% Pe = 120
for i=1:nz;
   ua120(i,1)=h5b(i+nz,2);
end;
%
% Pe = 180
for i=1:nz;
   ua180(i,1)=h5b(i+2*nz,2);
end;
%
% add series solution to plot
figure(2);
%
% Pe = 60
plot(z,u60(:,1),'-',z,ua60,'o');
hold on;
%
% Pe =120
plot(z,u120(:,1),'-',z,ua120,'+');
hold on;
%
% Pe = 180
plot(z,u180(:,1),'-',z,ua180,'*');
hold on;
%
% label axes
xlabel(' z');
ylabel(' u(0,z)');
title(' Graetz problem centerline temperature, solid - num, points - anal');
%
% print output
print h5p9b.eps -deps
%
% send Postscript file to ihd09 Postscript printer
!op -s 1 -d cs1 -q ihd109 h5p9a.eps
!op -s 1 -d cs1 -q ihd109 h5p9b.eps