function [t Y] = heun(F,Y0,T,h)
%Heun method for systems/or scalar equation
% input: 
%   F = @(t,X) [f_1(t,X(1),...,X(n)),...,f_n(t,X(1),...,X(n))]; function handle
%   Y0  row vector or a point, initial condition
%   T = [t0 tf]  row vector, t0 = starting time, tf = final time
%   h stepsize; default h = 0.02


%check input
if ~exist('h','var') || isempty(h)
    h = 0.02;
end

if (T(1) >= T(2))
    fprintf('Chyba: neplatny interval.\n');
    return
end

if h >= norm(T(2) - T(1),2)
    fprintf('Chyba: neplatne h.\n');
    return
end

%divide time interval
t = [T(1):h:T(2)]';
n = length(t);
if t(n) < T(2)
  t(n+1) = T(2);
elseif t(n) > T(2)
  t(n) = T(2);
  end
  
%initialize for approximate values
N = length(t);
vars = length(Y0);
Y = zeros(N,vars);


Y(1,:) = Y0;
m = 0;


for i = 1:(N-1)
    f1 = feval(F,t(i),Y(i,:));
    f2 = feval(F,t(i+1),Y(i,:) + h.*f1');
    m = m+ 2;
    Y(i+1,:) = Y(i,:) + h.*((f1' + f2')/2);%pozor na transpozici
end


[t Y];
fprintf('pocet vyhodnoceni F: %f\n',m)

  
end