# SECTION 3.7

# Example 2

# Remark: This solution is a modification of the solution to Example 1, Sect 3.6.

# We will use a Maple "procdure", using the proc command, to implement the subroutine on page 138. However, for simplicity, we won't implement the subroutine on page 139 (Runge-Kutta with tolerance).

f := (x,y) -> y;

NiM+JSJmR2orNiQlInhHJSJ5RzYiNiQlKW9wZXJhdG9yRyUmYXJyb3dHRik5JUYpRilGKTYkIiIhRi8=

x0 := 0; y0 := 1.; c := 1;

NiM+JSN4MEciIiE=

NiM+JSN5MEckIiIiIiIh

NiM+JSJjRyIiIg==

rk4 := proc(x0,y0,c,N,PRNTR)

local h,x,y,i,k1,k2,k3,k4:

h := evalf((c-x0)/N):

x := x0: y := y0:

for i from 1 to N do

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):

x := x+h:

y := y+1/6*(k1+2*(k2+k3)+k4):

if (PRNTR=1)

then print(x,y):

end if:

end do:

return(x,y):

end proc;

NiM+JSRyazRHais2JyUjeDBHJSN5MEclImNHJSJORyUmUFJOVFJHNiolImhHJSJ4RyUieUclImlHJSNrMUclI2syRyUjazNHJSNrNEc2IkY1Qyc+OCQtJSZldmFsZkc2IyomLCY5JiIiIjkkISIiRj85J0ZBPjglRkA+OCY5JT8oOCdGP0Y/RkIlJXRydWVHQyk+OCgqJkY4Rj8tJSJmRzYkRkRGRkY/PjgpKiZGOEY/LUZQNiQsJkZERj8qJiNGPyIiI0Y/RjhGP0Y/LCZGRkY/KiZGWUY/Rk1GP0Y/Rj8+OCoqJkY4Rj8tRlA2JEZXLCZGRkY/KiZGWUY/RlNGP0Y/Rj8+OCsqJkY4Rj8tRlA2JCwmRkRGP0Y4Rj8sJkZGRj9GaG5GP0Y/PkZERmNvPkZGLCxGRkY/KiYjRj8iIidGP0ZNRj9GPyomI0Y/IiIkRj9GU0Y/Rj8qJkZccEY/RmhuRj9GPyomRmlvRj9GX29GP0Y/QCQvOShGPy0lJnByaW50R0ZRT0ZRRjVGNUY1NiciIiFGZ3BGZ3BGZ3AiI2w=

# Before continuing, we check rk4 with the N=1 case:

rk4(x0,y0,c,1,0);

NiQkIiIiIiIhJCIrTExMM0YhIio=

# The above are the correct values of x and y.

# Next, we use a loop to change N:

Numbersteps := [1,2,4]:

print('N','y(1)');

NiQlIk5HLSUieUc2IyIiIg==

for k from 1 to 3 do

N := Numbersteps[k]:

out := rk4(x0,y0,c,N,0):

y := out[2]:

print(N,y);

end do:

NiQiIiIkIitMTEwzRiEiKg==

NiQiIiMkIishPll0ciMhIio=

NiQiIiUkIitRKjQjPUYhIio=

# Example 4

f := (x,y) -> x - y^2;

NiM+JSJmR2orNiQlInhHJSJ5RzYiNiQlKW9wZXJhdG9yRyUmYXJyb3dHRiksJjkkIiIiKiQpOSUiIiNGLyEiIkYpRilGKTYkIiIhRjY=

x0 := 0; y0 := 1.;c := 2;

NiM+JSN4MEciIiE=

NiM+JSN5MEckIiIiIiIh

NiM+JSJjRyIiIw==

Numbersteps := [1,2,4,8,16];

NiM+JSxOdW1iZXJzdGVwc0c3JyIiIiIiIyIiJSIiKSIjOw==

print('h','Approximation');

NiQlImhHJS5BcHByb3hpbWF0aW9uRw==

for k from 1 to 5 do

N := Numbersteps[k]:

y := rk4(x0,y0,c,N,0)[2]:

print(2./N,y);

end do:

NiQkIiIjIiIhJCErTUxMTCQpISIq

NiQkIisrKysrNSEiKiQiK3hrLnY3RiU=

NiQkIisrKysrXSEjNSQiK0BdcF43ISIq

NiQkIisrKysrRCEjNSQiKzotS143ISIq

NiQkIisrKytdNyEjNSQiK2ViSl43ISIq