fprint:=proc(ncase,neqn,t,u)
#
# Function fprint displays the numerical solution
# to the nonlinear PDE
#
# Declare global variables
  global nsteps:
#
# Type variables
  local
  im,             tmin,
  dx,              dxs,          i,
  alpha,           k,            sigma,
  L,               ui,           ua,
  a,               e,            rhos,
  cps,             ls,           cs:
#
# Problem parameters
  alpha:=1.0e-06: k:=1.0:        sigma:=5.67e-08:
  L:=0.1:         ui:=298.0:     ua:=2000.0:
  a:=1.0:         e:=1.0:        rhos:=7800.0:
  cps:=435.0:     ls:=0.025:     cs:=rhos*cps*ls:
#
# Print a heading for the solution at t = 0 
  if (t <= 0.0) then
#
#   Label for ODE integrator
#
#   Fixed step modified Euler 
    if (ncase = 1) then
      printf(`\n\n euler2a integrator\n\n`);
# 
#   Variable step modified Euler
    elif (ncase = 2) then
      printf(`\n\n euler2b integrator\n\n`);
#
#   Fixed step classical fourth order RK
    elif (ncase = 3) then
      printf(`\n\n rkc4a integrator\n\n`);
#
#   Variable step classical fourth order RK 
    elif (ncase = 4) then
      printf(`\n\n rkc4b integrator\n\n`);
#
#   Fixed step RK Fehlberg 45 
    elif (ncase = 5) then    
      printf(`\n\n rkc45a integrator\n\n`);
#        
#   Variable step RK Fehlberg 45 
    elif (ncase = 6) then
      printf(`\n\n rkc45b integrator\n\n`);
    end if:
#
# Heading 
  printf(` ncase = %2d   neqn = %2d   nsteps = %3d\n\n`,ncase,neqn,nsteps);
  printf(`    t      u(1)      u(im)     u(n-1)\n`);        
#
# End of t = 0 heading
  end if:
#
# Numerical solution output
#   
# Grid index of midpoint 
  im:=(neqn/2);
# 
# Display the numerical solution
  tmin:=t*L^2/alpha/60.0;
  printf(`%7.1f%10.2f%10.2f%10.2f\n`,tmin,u[1],u[im],u[neqn-1]);
#
# End of fprint
  end: