function [X,Y]=taylorova_metoda(f, y0, x0, xn, h, RAD, 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.
% Metoda ma rad = RAD  

is_symbolic=is_toolbox('symbolic');

if ~is_symbolic
  clc;
  disp(' ');
  disp('Lituji, ale neni k dispozici Symbolic Math Toolbox !');
  disp(' '); 
else

if zobraz==1
 clc;
 disp('----------------------------------------------------------------');
 disp(sprintf('  METODA TAYLOROVA TYPU radu %d pro reseni',RAD));
 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.  ');
end;

if RAD>1
  syms x y;
  derivacef=diff(f,x)+diff(f,y)*f;
  F=[f;derivacef];
  for j=2:RAD-1
    der=diff(F(j),x)+diff(F(j),y)*f;
    F=[F;der];
  end;
end;


n=ceil((xn-x0)/h);
X=x0:h:x0+n*h;
Y=zeros(size(X));
Y(1)=y0;
if RAD==1
  for i=1:n	
    x=x0+(i-1)*h;
    y=Y(i);
    hodnota=eval(f);
    Y(i+1)=y+h*hodnota;   
  end;
else
  for i=1:n	
    x=x0+(i-1)*h;
    y=Y(i);
    Y(i+1)=y;
    for j=1:RAD
      hodnota=eval(F(j));
      Y(i+1)=Y(i+1)+h^j/factorial(j)*hodnota;
    end;
  end;
end;

if zobraz==1
 disp(' ');
 disp('  Stiskni libovolnou klavesu.');
 disp(' ');
 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               ');
 popis=sprintf('Metoda Taylorova typu radu %d',RAD);
 title(popis);
 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;

end;







