Euler := proc( r, y0, n, f) local x0, x1, x, y, h, i; # internal variables x0 := op(1,r); x1 := op(2,r); # endpoints of the interval h := evalf((x1 - x0)/n); # compute step size x := x0; y := y0; # initialize for i from 1 to n do y := EulerStep( x, y, h, f ); # compute new y-value x := x + h; # compute new x-value od; RETURN(y); # value of Euler(...) end: ========================================================== Eulerplot := proc( r, y0, n, f) local x0, x1, i, x, y, h, data; # internal variables x0 := op(1,r); x1 := op(2,r); # endpoints of interval h := evalf((x1 - x0)/n); # compute step size x := x0; y := y0; # initial point data := array(0..n); # array to hold the data data[0] := [x0,y0]; for i from 1 to n do y := EulerStep( x, y, h, f ): # compute new y-value x := x + h: # compute new x-value data[i] := [x,y] od: plot( convert( data, list ) ); # convert to list & plot end: ========================================================== ndsolve := proc( r, y0, n, f, G) local x0, x1, x, y, h, i; # internal variables x0 := op(1,r); x1 := op(2,r); # endpoints h := evalf((x1 - x0)/n); # step size x := x0; y := y0; # initialize for i from 1 to n do y := G( x, y, h, f ); # new y-value x := x + h; # new x-value od; RETURN(y); # value of Euler(...) end: ========================================================== ndeplot := proc( r, y0, n, f, G) local x0, x1, i, x, y, h, data; # internal variables x0 := op(1,r); x1 := op(2,r); # endpoints of interval h := evalf((x1 - x0)/n); # compute step size x := x0; y := y0; # initial point data := array(0..n); # array to hold data data[0] := [x0,y0]; for i from 1 to n do y := G( x, y, h, f ): # compute new y-value x := x + h: # compute new x-value data[i] := [x,y] od: plot( convert( data, list ) ); # convert to list & plot end: ========================================================== EulerStep := ( x0, y0, h, f ) -> y0 + h*f( x0,y0 ); ========================================================== HeunStep := proc( x, y, h, f) local m; m := f( x, y ); y+h*( m + f( (x+h), (y+h*m) ))/2; end: ========================================================== rkStep := proc(x, y, h, f) local k1, k2, k3, k4; k1 := h*f(x, y); k2 := h*f( (x + h/2), (y + k1/2) ); k3 := h*f( (x + h/2), (y + k2/2) ); k4 := h*f((x +h), (y + k3) ); y + ( k1 + 2*k2 + 2*k3 + k4)/6; end: ==========================================================