Mathematical Principles for Scientific Computing and Visualization
http : // www.farinhansford.com/books/scv
by Gerald Farin and Dianne Hansford

This notebook demonstrates Heun' s method for solving an initial value ODE problem.
Problem: Find values of y over [0,2] for y' = y given initial value y(0) = 1.
The exact solution is e^x.

Heun's method is allowed 10 steps. Compare this to the solution using Euler's method in Euler.nb

The figure generated here and the method are described in detail in Chapter 9 : Computing Dynamic Processes

Updated May 2008

(* Set output foleder to be folder this files lives in *)
(SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];) 


(* Huen's method *)
heun[f_,{x_,x0_,xn_},{y_,y0_},steps_]:=Block[{xold=x0,yold=y0,sollist={{x0,y0}},x,y,h},h=N[((xn-x0)/steps)];Do[xnew=xold+h;ynew1=yold+h (f/.{x->xold,y->yold});oldslope=f/.{x->xold,y->yold};newslope=f/.{x->xnew,y->ynew1};ynew=yold+(1/2) h (oldslope+newslope);sollist=Append[sollist,{xnew,ynew}];xold=xnew;yold=ynew,{steps}];Return[sollist]]


(* Store the plot the exact solution e^x *)
exactplot=Plot[(E^x),{x,0,2},PlotStyle->{GrayLevel[0.7],Thickness[0.03`]},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}, AxesStyle->Directive[GrayLevel[0.0],Thickness[0.005]],TicksStyle->Directive[GrayLevel[0.0],Thickness[0.005]], ColorOutput->CMYKColor];


(* Solve with Huen's for 10 steps *)
heunsol=heun[y,{x,0,2},{y,1.`},10];

(* Store the plot Huen's solution *)
heunplot=ListPlot[heunsol,PlotStyle->{GrayLevel[0],PointSize[0.03`]},Frame->True,BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}, AxesStyle->Directive[GrayLevel[0.0],Thickness[0.005]],TicksStyle->Directive[GrayLevel[0.0],Thickness[0.005]], ColorOutput->CMYKColor];

(* Plot exact and Huen's together *) 
plt=Show[{exactplot,heunplot}]

(* Send to postscript file *)
Export["heun10.eps",plt]