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('----------------------------------------------------------------');
 fprintf('  METODA TAYLOROVA TYPU radu %d pro reseni\n',RAD);
 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.  ');
end

if RAD>1
  x = sym('x');
  y = sym('y');
  fev=eval(f);
  derivacef=diff(fev,x)+diff(fev,y)*fev;
  FFF=[fev;derivacef];
  for j=2:RAD-1
    der=diff(FFF(j),x)+diff(FFF(j),y)*fev;
    FFF=[FFF;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	
    xx=x0+(i-1)*h;
    yy=Y(i);
    Y(i+1)=yy;
    for j=1:RAD
      %hodnota=eval(FFF(j));
      hodnota=double(subs(FFF(j),[x,y],[xx,yy]));
      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,:),'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

end







