function d = dum()

% SAVED FROM THE DEATH

% Convergentie in functie van h
%
% voor het oplossen van volgend beginwaardeprobleem
%
%     y' = f(t,y)    met    y(a) = y0
%     f(t,y) = -2 t y^2 in [a,b] = [0,5]   y0 = 1
%
% fout in y(b) in functie van de stap h
%
% Opgelet, de files f.m en fexact.m zullen worden gecreerd en de bestaande
% files met deze naam zullen worden overschreven!

clf; clear; format long; echo on;

if(exist('f.m') == 2), delete f.m; end;
diary f.m; disp('function z = f(t,y);');...
           disp('% het rechterlid Dy= f(t,y)');...
           disp('z = -2*t*y*y;');...
diary off; disp(' '); exist('f.m');

if(exist('fexact.m') == 2), delete fexact.m; end;
diary fexact.m; disp('function z = fexact(t);');...
                disp('% de exacte oplossing');...
                disp('z = 1./(1+t.^2);');...
diary off; disp(' '); exist('fexact.m');

disp('Druk een willekeurige toets om verder te gaan.'); pause

%a  = input('Geef beginpunt. ');
%b  = input('Geef eindpunt. ');
%ya = input('Geef beginwaarde y(a). ');
% m  = input('Geef het aantal deelintervallen van [a,b]. ');
a=0; b=5; ya=1;

% exacte oplossing
exac=feval('fexact',b);
col='rbg';
hold off;
for meth = 1:3,
 EE=[]; TT=[];
 for m= 1:10,
  if meth == 1
	 [T,Y] = feuler('f',a,b,ya,10*m);
  elseif meth == 2
	 [T,Y] = fheun('f',a,b,ya,10*m);
  else
	 [T,Y] = fcauchy('f',a,b,ya,10*m);
  end
  i=size(Y);
  EE=[EE abs(Y(i(2))-exac)]; TT=[TT (b-a)/(10*m)];
  errpl=EE;
 end
 if max(errpl) ~= 0
   ipl=find(errpl == zeros(size(errpl)));
   if length(ipl) == 1
     errpl(ipl)=eps;
   else
     errpl(ipl)=eps*ones(size(ipl));
   end
 end
 semilogy(errpl,col(meth)), title('absolute fout in functie van de stap')
 hold on
end
% print -depsc cvg.eps

function [T,Y] = feuler(f,a,b,ya,m)
% [T,Y] = feuler(f,a,b,ya,m)
% Forward Euler's solution for y' = f(t,y) with y(a) = ya.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% ya is the initial condition, input.
% m  is the number of steps, input.
% T  is the vector of abscissas, output.
% Y  is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
  Y(j+1) = Y(j) + h*feval(f,T(j),Y(j));
  T(j+1) = a + h*j;
end


function [T,Y] = fcauchy(f,a,b,ya,m)
% [T,Y] = fcauchy(f,a,b,ya,m)
% Cauchy's solution for y' = f(t,y) with y(a) = ya.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% ya is the initial condition, input.
% m  is the number of steps, input.
% T  is the vector of abscissas, output.
% Y  is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
  k1=feval(f,T(j),Y(j));
  k2=feval(f,T(j)+h/2,Y(j)+h*k1/2);
  Y(j+1) = Y(j) + h*k2;
  T(j+1) = a + h*j;
end

function [T,Y] = fheun(f,a,b,ya,m)
% [T,Y] = feuler(f,a,b,ya,m)
% Forward Euler's solution for y' = f(t,y) with y(a) = ya.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% ya is the initial condition, input.
% m  is the number of steps, input.
% T  is the vector of abscissas, output.
% Y  is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
  k1=feval(f,T(j),Y(j));
  k2=feval(f,T(j)+h,Y(j)+h*k1);
  Y(j+1) = Y(j) + 0.5*h*(k1+k2);
  T(j+1) = a + h*j;
end