#############################################################
# Building 'bvp1d' and sample output on the SGIs
#############################################################

einstein% ls Makefile  bvp1d.f
Makefile   bvp1d.f

einstein% make
        f77 -g -n32 -c bvp1d.f
        f77 -g -n32 -L/usr/localn32/lib -n32 bvp1d.o \
            -lp329f -llapack -lblas -o bvp1d

einstein% bvp1d
 usage: bvp1d <level> [<option>]

        Specify option .ne. 0 for output
        of error instead of solution

einstein% bvp1d 4
   0.0000000000000000E+00 -5.5511151231257827E-16
   6.2500000000000000E-02  0.3876394685723090    
   0.1250000000000000      0.7162643420150174    
   0.1875000000000000      0.9358444623383684    
   0.2500000000000000       1.012950746721879    
   0.3125000000000000      0.9358444623383684    
   0.3750000000000000      0.7162643420150175    
   0.4375000000000000      0.3876394685723092    
   0.5000000000000000     -2.2204460492503131E-16
   0.5625000000000000     -0.3876394685723097    
   0.6250000000000000     -0.7162643420150181    
   0.6875000000000000     -0.9358444623383690    
   0.7500000000000000      -1.012950746721880    
   0.8125000000000000     -0.9358444623383690    
   0.8750000000000000     -0.7162643420150181    
   0.9375000000000000     -0.3876394685723097    
    1.000000000000000     -2.4492935982947064E-16
 rmserr =   8.8841389573651453E-03

#############################################################
# Simple convergence test: solve BVP on a sequence of 
# levels (h, h/2, h/4, h/16, etc.), redirect stdout to 
# /dev/null so that only the overall RMS error appears on 
# terminal.  Note how RMS error goes down by very nearly 
# a factor of 4 at each successive level, indicating 
# O(h^2) convergence.
#############################################################
einstein% foreach level (4 5 6 7 8 9 10)
foreach? bvp1d $level > /dev/null
foreach? end
 rmserr =   8.8841389573651453E-03
 rmserr =   2.2413991373367772E-03
 rmserr =   5.6382739826354859E-04
 rmserr =   1.4145099550532311E-04
 rmserr =   3.5428279660444339E-05
 rmserr =   8.8654982501522291E-06
 rmserr =   2.2174426911240527E-06

#############################################################
# Making output files for subsequent plotting via gnuplot.
# See Class Notes for postscript.
#############################################################
einstein% bvp1d 8 > out8
 rmserr =   3.5428279660453927E-05

einstein% bvp1d 5 1 > err5
 rmserr =   2.2413991373375825E-03

einstein% bvp1d 6 1 > err6
 rmserr =   5.6382739826356723E-04

einstein% bvp1d 7 1 > err7
 rmserr =   1.4145099550242024E-04

#############################################################
# Gnuplot 'script' (gpsoln8) for making plot of level-8 
# solution
#############################################################
einstein% cat gpsoln8

set terminal postscript portrait
set title "Solution of 2-pt BVP using O(h^2) FDA   (level 8, 257 mesh points)"
set xlabel "x"
set ylabel "u(x)"
set size 0.760,1.0
set output "soln8.ps"
plot [0:1] [-1:1] "out8" t "u(x)"
quit

#############################################################
# Make the plot
#############################################################
einstein% gnuplot < gpsoln8

#############################################################
# Gnuplot 'script' (gperr) for making plot of error from 
# level 5, 6 and 7 calculations
#############################################################
einstein% cat gperr

set terminal postscript portrait
set title "Error in solution of 2-pt BVP using O(h^2) FDA  (levels 5, 6, 7)"
set xlabel "x"
set ylabel "u_exact(x) - u(x)"
set size 0.760,1.0
set output "err567.ps"
plot "err5" t "level 5", "err6" t "level 6", "err7" t "level 7"
quit

#############################################################
# Make the plot
############################################################
einstein% gnuplot < gperr