function [X,Y]=heunova_metoda(f, y0, x0, xn, h, zobraz);
%
% Vypocte reseni obycejne diferencialni rovnice 1. radu
%                      y'=f(x,y(x))
% s pocatecni podminkou
%                       y(x0)=y0
% na intervalu  <x0,xn>.
% Pro diskretizaci pouzije krok h 

if zobraz==1
 clc;
 disp('----------------------------------------------------------------');
 disp(sprintf('  HEUNOVA METODA pro reseni'));
 disp(' ');
 disp(sprintf('  obycejne diferencialni rovnice 1. radu'));
 disp(sprintf('                    y`=%s',f));
 disp(sprintf('  s pocatecni podminkou'));
 disp(sprintf('                    y(%f)=%f',x0,y0));
 disp(sprintf('  na intervalu'));
 disp(sprintf('                    < %f, %f >',x0,xn));
 disp(sprintf('  s krokem h=%f',h));
 disp('----------------------------------------------------------------');
 disp(' ');
 disp('  Chvilku strpeni, probiha vypocet.  ');
 disp(' ');
end;

n=ceil((xn-x0)/h);
X=x0:h:x0+n*h;
Y=zeros(size(X));
Y(1)=y0;
for i=1:n	
  x=x0+(i-1)*h;
  y=Y(i);
  k1=eval(f);
  x=x+h;
  y=y+h*k1;
  k2=eval(f);
  Y(i+1)=Y(i)+h*(k1+k2)/2;
end;

if zobraz==1
 disp(' ');
 disp('  Stiskni libovolnou klavesu.');
 max_y=max(Y);
 min_y=min(Y);
 sy=max_y-min_y;
 sx=X(length(X))-X(1);
 figure(1);
 clf;
 hold on;
 zoom on;
 axis([X(1)-.15*sx X(length(X))+.15*sx min_y-.15*sy max_y+.15*sy]);
 texty(1,:)=sprintf(' Pocatecni podminka  ');
 texty(2,:)=sprintf(' Priblizne reseni    ');
 title('Heunova metoda');
 plot(X(1),Y(1),'go');
 legend(texty(1,:),0);
 pause;
 for i=2:n+1
   plot(X(i-1:i),Y(i-1:i),'b-');
   if i==2
     legend(texty,0);
   end;
   plot(X(i),Y(i),'ro');
   if i<n+1
     pause
   end;
 end;
end;









