function [t,y]=linear8(fcn,t0,tend,y0,y1,steps);
%-----------------------------------------------
% INPUT
% fcn   : function describing the IVP
% t0    : left side of integration
% tend  : right side of integration
% y0    : column vector of y(t0)
% y1    : column vector of y(t0+h)
% steps : steps of integration, h=(tend-t0)/steps

% OUTPUT
% t   : array containing the values of t
% y   : array of size n X steps with the values of y
%-----------------------------------------------

c1=-0.4522764562357594;c2=0.454777110594506;
c3=-0.8230131465004396;c4=0.83205029505874633;
d32=0.6044425455247441;w1=-7.176517980971427e-004;
w2=0.420891641969056;b1=0.238872506910894;
b2=0.240505260009359;b3=0.050874312458376;
b4=0.049573930450411;d11=0.388829320534869;
d12=-0.060414095984401;d21=-0.378268263774167;
d22=0.228651575950437;a21=0.0256392426868207;
d31=-1.495977126547815;d32=0.604442545524744;
a31=-1.142558822542608;a32=2.784275296472176;
d41=1.497296269018007;d42=-0.296348071675540;
a41=1.316152612845400;a42=-2.588722480460740;
a43=0.001750369497174;

y=zeros(length(y0),steps+1);t=zeros(steps+1,1);
h=(tend-t0)/steps;x=h;
t(1)=t0;y(:,1)=y0;t(2)=t0+h;y(:,2)=y1;
f0=feval(fcn,t0,y0);f1=feval(fcn,t0+h,y1);

for i1=2:steps,
ya1=(d11*f0+d12*f1)*h^2+c1*y0+(1-c1)*y1;
fa1=feval(fcn,x-c1*h,ya1);

ya2=(d21*f0+d22*f1+a21*fa1)*h^2+c2*y0+(1-c2)*y1;
fa2=feval(fcn,x-c2*h,ya2);

ya3=(d31*f0+d32*f1+a31*fa1+a32*fa2)*h^2+c3*y0+(1-c3)*y1;
fa3=feval(fcn,x-c3*h,ya3);

ya4=(d41*f0+d42*f1+a41*fa1+a42*fa2+a43*fa3)*h^2+c4*y0+(1-c4)*y1;
fa4=feval(fcn,x-c4*h,ya4);

y2=(w1*f0+w2*f1+b1*fa1+b2*fa2+b3*fa3+b4*fa4)*h^2-y0+2*y1;
x=x+h;y0=y1;y1=y2;f0=f1;f1=fcn(x,y1);
t(i1+1)=x;y(:,i1+1)=y2;
end


% Example from the paper
% the function of dimension n=2
% fcn=@(t,y) [.01*y(1)-.1*y(2);-.1*y(1)+.01*y(2)+sin(t)];

% initial values
% y0=[1 1]';
% h=pi/24;
% y1=[cos(.3*h)-1000/10101*sin(h) cos(.3*h)-10100/10101*sin(h)]';

% the integration
% [t,y]=linear8(fcn,0,10*pi,y0,y1,240);

% end point error (observe that y(10*pi)==y(:,241)=[-1 -1]')
% max(-log10(abs(y(:,241)+[1 1]')))
% ans =
%   10.9786