Dear all;

I want to plot the error defined in the last line. the error is defined by the difference entre step h and step 2h for Heun method.

The code for Heun1 work well ( stepsize 2h),

Heun2 don't working. Please, I need a Help, maybe there is somethong missing.

here the code with maple code, without pciture.

f :=(x,y)-> 1/(1+sin(y));

ode:=diff(y(x),x)=f(x,y);

## Heun Method with stepsize 2*h

Heun1 := proc(f, x0, h)

local x, y, i, N, k;

N := round((1/2)*x0/h);

y := Array(0 .. N);

x := Array(0 .. N);

x[0] := 0;

y[0] := (1/4)*Pi;

for i from 0 to N-1 do

x[i+1] := (i+1)*2*h;

k[1] := f(x[i], y[i]);

k[2] := f(x[i]+h, y[i]+h*k[1]);

y[i+1] := y[i]+(1/2)*h*(k[1]+k[2]);

end do;

evalf([seq([x[i], y[i]], i = 0..N)])

end proc

## Heun Method with stepsize h

Heun2 := proc(f, x0, h)

local x, y, i, N, k;

N := round((1/2)*x0/h);

y := Array(0 .. 2*N);

x := Array(0 .. 2*N);

x[0] := 0;

y[0] := (1/4)*Pi;

for i from 0 to 2*N-1 do

x[i+1] := (i+1)*h;

k[1] := f(x[i], y[i]);

k[2] := f(x[i]+h, y[i]+h*k[1]);

y[i+1] := y[i]+(1/2)*h*(k[1]+k[2]);

end do;

evalf([seq([x[2*i], y[2*i]], i = 0..N)])

end proc

Heun1((x,y)->f(x, y), 1,1e-2); ## Heun1 work very well

Heun2((x,y)->f(x, y), 1,1e-2); ## Dont working...

## The error

x0:=1;

epsilon:= (x0,h)->sqrt(1/(( round(x0/(2*h))+1))*add(( Heun1((x,y)->f(x,y),x0,round(x0/(2*h)))[i][2]-Heun2((x,y)->f(x,y),x0, round(x0/(2*h)))[i][2])^(2),i=1..N+1));

data[Heun] := [seq([h,epsilon(x0,h)], h =1e-2..1)];

loglog(epsilon(1,h),h=1e-2..1);