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