Hi!

I want to write a Mathematica program for Heun's Method. I am trying to modify some code I found of Euler's method:

euler[f_, {t_, t0_, tn_}, {y_, y0_}, steps_] :=

Block[{told = t0, yold = y0, sollist = {{t0, y0}}, t, y, h},

h = N[(tn - t0)/steps];

Do[

tnew = told + h;

ynew = yold + h*(f /. {t -> told, y -> yold});

sollist = Append[sollist, {tnew, ynew}

];

told = tnew;

yold = ynew, {steps}];

Return[sollist]]

The Line I think needs to be modified is:

ynew = yold + h*(f /. {t -> told, y -> yold});

I came up with:

ynew = yold + (h/2)*(f /. {t -> told, y -> yold} + (f /. {t -> told + h, y -> yold + h*f /. {t -> told, y -> yold}} ));

Which usually makes either Mathematica or my computer crash.

I could really use some help, I'm not very experienced with Mathematica.

I hope I have provided enough info, and thanks in advance!

Links to helpful wiki articles:

Euler method - Wikipedia, the free encyclopedia

Heun's method - Wikipedia, the free encyclopedia