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('----------------------------------------------------------------');
 fprintf('  HEUNOVA METODA pro reseni\n');
 disp(' ');
 fprintf('  obycejne diferencialni rovnice 1. radu\n');
 fprintf('                    y`=%s\n',f);
 fprintf('  s pocatecni podminkou\n');
 fprintf('                    y(%f)=%f\n',x0,y0);
 fprintf('  na intervalu\n');
 fprintf('                    < %f, %f >\n',x0,xn);
 fprintf('  s krokem h=%f\n',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,:),'Location','bestoutside');
 pause;
 i=2;
   plot(X(i-1:i),Y(i-1:i),'b-');
   legend(texty,'Location','bestoutside');
   bbb=plot(X(i),Y(i),'ro');
   set(get(get(bbb,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
   pause
 for i=3:n+1
   bbb=plot(X(i-1:i),Y(i-1:i),'b-');
   set(get(get(bbb,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
   bbb=plot(X(i),Y(i),'ro');
   set(get(get(bbb,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
   if i<n+1
     pause
   end
 end
end









