function [yn1,en]=rkf(f,tn,hn,yn,a,b,c,est)
%
%     Copyright: Eduardo Casas (Octubre, 2000)
%
m=length(b);
kn=feval(f,tn,yn);
for i=2:m
   kn=[kn,feval(f,tn+hn*c(i),yn+hn*(kn*a(i,1:i-1)'))];
end
yn1=yn+hn*(kn*b);
en=hn*norm(kn*est);
