% problem41
function problem41
resul=zeros(5,2);
irep=1;
for tol=[1e-6,1e-7,1e-8,1e-9],
    [tout, yout, dyout, fev] = grkn(@(x,y,dy) -5*dy(1)-y(1)+sin(x/10), 0,10, 0, 0, tol);
     resul(irep,:)=[fev max(abs(yout(:,end)-0.50814725856006851284))];
     irep=irep+1;
end;

axis([150 900 4 12]);
loglog(resul(:,1),resul(:,2),'-.^k');
xlabel('function evaluations');
ylabel('accurate digits');
legend('NEW7(5)');
return


% various grkn 
function [tout, yout, dyout, fev] = grkn(FunFcn, t0, tfinal, y0, dy0, tol);

%-------------------------------------------------------
% my new method 7(5)
c=[0,1/8,1/5,2/5,1/2,3/5,4/5,5/6,1]';
b=[0.20119597167401398108,-1.2308182696272720194,2.2061162664677621851,-2.8566625938490179854,3.3864387812136187983,-0.94554726562324883883,-0.90016335789451293256,1.1394404676386568117,0];
bb=[-0.23104875124991469820,1.5381600865012839022,-1.3057711936930409504,0.15338322099337164246,0.80266388293964741919,-0.22271536117021783691,-0.31164355085541809419,0.52697166653428861587,0.05];
a=[[0,0,0,0,0,0,0,0,0],
[0.125,0,0,0,0,0,0,0,0],
[0.088447245894008380024,0.11155275410599161998,0,0,0,0,0,0,0],
[-0.0057453039845693144462,-0.23251345088521661222,0.63825875486978592666,0,0,0,0,0,0],
[-0.057122827093073277786,-0.25416461483902363251,0.67753596559208772069,0.13375147634000918960,0,0,0,0,0],
[0.14575738023521525598,-1.1046796859421592459,1.5972258006178789110,-0.43100128372938146827,0.39269778881844654711,0,0,0,0],
[0.22855740566333236569,0.60553735393465262555,-0.55119449835624047120,0.37446628565273008921,-0.63306953205127876816,0.77570298515680415892,0,0,0],
[0.44204418043038272803,-0.27328844564647506815,0.42387844189337618372,-0.28790724390541638102,-0.14873885875264402674,0.62398772886280406679,0.053357530451305830706,0,0],
b];
ahat=[[0, 0, 0, 0, 0, 0, 0, 0, 0], 
      [0, 0, 0, 0, 0, 0, 0, 0, 0], 
      [0.01394409426324895250, 0, 0, 0, 0, 0, 0, 0, 0], 
      [0.02738804767531949790, 0.07119952193798561090, 0, 0, 0, 0, 0, 0, 0], 
      [0.02738717040592154645, 0.04448198564285182197, 0.08536805075076989749, 0, 0, 0, 0, 0, 0], 
      [-0.01677051210500938968,0.1785786505607719841, -0.009023467167410410290, 0.05252390900992437864, 0, 0, 0, 0, 0], 
      [0.1740162676415949263, -0.8445551689042742859, 1.049051830044885967, -0.4190029669330332324, 0.3046168470509452561, 0, 0, 0, 0], 
      [0.1166263992101645447, -0.5049653213888277587, 0.6827036772788881907, -0.2688529578788410725, 0.2112595745400474752, 0.04138959565167301271, 0, 0, 0], 
      [0.02436559579726689012, 0.2376463188860853530, -0.05995898849597635318, 0.1953392780892177322, 0.02907260054525608803, 0.01273746572867050314, 0.06079772944947978678, 0, 0]];
bhat=[0.024365595797266890,0.237646318886085353,-0.05995898849597635,0.195339278089217732,0.029072600545256088,0.012737465728670503,0.06079772944947978678,0,0];
bbhat=[0.169360189502373042,-0.53357364861436657,0.79145906450316170,-0.207904006838439304,0.2007729818247841318,0.0398036573912011186,-0.01689026115064696005,0.05697202338193284059,0];
pow=1/6;
%-------------------------------------------------------

if nargin < 5, tol = 1.e-6; end

% Initialization
t = t0;
hmax = (tfinal - t)/5;
hmin = (tfinal - t)/2000000;
y = y0(:);
dy = dy0(:);
f = y*zeros(1,length(c));
tout = t;
yout = y;
dyout=dy;
fev = 1;

f(:,1) = feval(FunFcn,t,y,dy);
h=tol^pow/max(max(abs(f(:,1))),1);

% The main loop
   while (t < tfinal) && (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes
      f(:,1) = feval(FunFcn,t,y,dy);
      for j = 2:length(c),
         f(:,j) = feval(FunFcn, t+c(j)*h, y+c(j)*h*dy+h^2*(ahat(j,:)*f')', dy+h*(a(j,:)*f')');
      end



      % Estimate the error and the acceptable error
      delta = max(norm(h^2*(bhat-bbhat)*f','inf'), norm(h*(b-bb)*f','inf'));

      % Update the solution only if the error is acceptable
      if delta <= tol,
         t = t + h;
         y = y + h*dy + h^2*(bhat*f')';
         dy = dy + h *(b*f')';
         tout = [tout t];
         yout = [yout y];
         dyout = [dyout dy];
         fev = fev+length(c)-1;
      else
         fev = fev+length(c)-1; 
      end

      % Update the step size
      if delta ~= 0.0
         h = min(hmax, 0.9*h*(tol/delta)^pow);
      end
   end;

   if (t < tfinal)
      disp('SINGULARITY LIKELY.')
      t
   end