trapetsi := proc (t0,tf,n,x0) 
	local a, b , tv,xv, h, ap0,ap1,k,xk,tk; 
    h := evalf((tf-t0)/n); 
	tv := Vector(n+1); xv := Vector(n+1);
	tv[1] := t0; xv[1] := x0; 
	for k to n do 
		xk:=xv[k];tk:=tv[k];
		tv[k+1]:=tk+h;
		ap0:=xk^2-tk^2-(tk+h)^2;
		ap1:=xk+ap0*h/2;
		xv[k+1] := (1-sqrt(1-2*h*ap1))/h;
	end do ;
	Matrix([tv,xv]);
end proc;



eulerinmen := proc (fu,t0,tf,n,x0) 
	local a, b , tv,xv, h,k,xk,tk; 
    h := evalf((tf-t0)/n); 
	tv := Vector(n+1); xv := Vector(n+1);
	tv[1] := t0; xv[1] := x0; 
	for k to n do 
		xk:=xv[k];tk:=tv[k];
		tv[k+1]:=tk+h;
		xv[k+1] := xk+h*fu(tk,xk);
	end do ;
	Matrix([tv,xv]);
end proc;





heuninmen := proc (fu,t0,tf,n,x0) 
	local a, b , tv,xv, h,k,xk,tk,k1,k2; 
    h := evalf((tf-t0)/n); 
	tv := Vector(n+1); xv := Vector(n+1);
	tv[1] := t0; xv[1] := x0; 
	for k to n do 
		xk:=xv[k];tk:=tv[k];
		k1:=fu(tk,xk);
		k2:=fu(tk+h,xk+h*k1);
		tv[k+1]:=tk+h;
		xv[k+1] := xk+h*(k1+k2)/2;
	end do ;
	Matrix([tv,xv]);
end proc;


ruku4 := proc (fu,t0,tf,n,x0) 
	local a, b , tv,xv, h,k,xk,tk,k1,k2,k3,k4; 
    h := evalf((tf-t0)/n); 
	tv := Vector(n+1); xv := Vector(n+1);
	tv[1] := t0; xv[1] := x0; 
	for k to n do 
		xk:=xv[k];tk:=tv[k];
		k1:=fu(tk,xk);
		k2:=fu(tk+h/2,xk+h*k1/2);
		k3:=fu(tk+h/2,xk+h*k2/2);
		k4:=fu(tk+h,xk+h*k3);
		tv[k+1]:=tk+h;
		xv[k+1] := xk+h*(k1+2*k2+2*k3+k4)/6;
	end do ;
	Matrix([tv,xv]);
end proc;





virhe:=proc(fu,ma)
	local m,vi,t,x,j ;
	m:=RowDimension(ma);
	vi:=Matrix(m,2);
	for j to m do
		t:=ma[j,1];
		x:=ma[j,2];
		vi[j,1]:=t;
		vi[j,2]:=fu(t)-x;
	end do;
	vi;
end proc;

virhemax:=proc(fu,ma)
	local m,vi,t,x,j ;
	m:=RowDimension(ma);
	vi:=Vector(m);
	for j to m do
		t:=ma[j,1];
		x:=ma[j,2];
		vi[j]:=fu(t)-x;
	end do;
	max(abs(vi));
end proc;


eusup:=proc (fu,t0,tf,n,x0,m,f)
	local va,su,j ,ap,k;
	su:=Matrix(m+1,2);
	for j from 0 to m do
		k:=n*2^j;
		va:=eulerinmen(fu,t0,tf,k,x0);
		ap:=virhemax(f,va);
		su[j+1,1]:=ap;
		if (j>0) then
			su[j,2]:=su[j,1]/su[j+1,1];
		end if;
	end do;
	su;
end proc;


eusup1:=proc (fu,t0,tf,n,x0,m,f)
	local va,su,j ,ap,k,h;
	su:=Matrix(m+1,2);
	for j from 0 to m do
		k:=n*2^j;
		h:=(tf-t0)/k;
		va:=eulerinmen(fu,t0,tf,k,x0);
		ap:=virhemax(f,va);
		su[j+1,1]:=h;
		su[j+1,2]:=ap;
	end do;
	pointplot(su, color = "DarkGreen", symbolsize = 20, symbol = solidcircle, axis = [mode = log]);
end proc;

hesup:=proc (fu,t0,tf,n,x0,m,f)
	local va,su,j ,ap,k;
	su:=Matrix(m+1,2);
	for j from 0 to m do
		k:=n*2^j;
		va:=heuninmen(fu,t0,tf,k,x0);
		ap:=virhemax(f,va);
		su[j+1,1]:=ap;
		if (j>0) then
			su[j,2]:=su[j,1]/su[j+1,1];
		end if;
	end do;
	su;
end proc;

hesup1:=proc (fu,t0,tf,n,x0,m,f)
	local va,su,j ,ap,k,h;
	su:=Matrix(m+1,2);
	for j from 0 to m do
		k:=n*2^j;
		h:=(tf-t0)/k;
		va:=heuninmen(fu,t0,tf,k,x0);
		ap:=virhemax(f,va);
		su[j+1,1]:=h;
		su[j+1,2]:=ap;
	end do;
	pointplot(su, color = "Red", symbolsize = 20, symbol = solidcircle, axis = [mode = log]);
end proc;


rk4sup:=proc (fu,t0,tf,n,x0,m,f)
	local va,su,j ,ap,k;
	su:=Matrix(m+1,2);
	for j from 0 to m do
		k:=n*2^j;
		va:=ruku4(fu,t0,tf,k,x0);
		ap:=virhemax(f,va);
		su[j+1,1]:=ap;
		if (j>0) then
			su[j,2]:=su[j,1]/su[j+1,1];
		end if;
	end do;
	su;
end proc;


rk4sup1:=proc (fu,t0,tf,n,x0,m,f)
	local va,su,j ,ap,k,h;
	su:=Matrix(m+1,2);
	for j from 0 to m do
		k:=n*2^j;
		h:=(tf-t0)/k;
		va:=ruku4(fu,t0,tf,k,x0);
		ap:=virhemax(f,va);
		su[j+1,1]:=h;
		su[j+1,2]:=ap;
	end do;
	pointplot(su, color = "DarkGoldenrod", symbolsize = 20, symbol = solidcircle, axis = [mode = log]);
end proc;