function [t,y] = RDE43(fode,tspan,y0,options);

%   Reserch code
%
%   RDE43 solves stiff ordinary differential equations arising from parabolic
%   problems by Exponential W-methods [1,3]. The matrix functions inside the method
%   are computed using the Restricted-Denominator Rational Arnoldi method [2].
%
%   [T,Y] = RDE43(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the
%   system of differential equations y' = f(y) from time T0 to TFINAL with
%   initial conditions Y0. Function ODEFUN(T,Y) must return a column vector
%   corresponding to f(y). Each row in the solution array Y corresponds to
%   a time returned in the column vector T.
%
%   [T,Y] = RDE43(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
%   integration properties replaced by values in OPTIONS. RDE43 supports the
%   most important features of the matlab ODE solver. Moreover in this code
%   there is another option available [3]:
%
%       'Method'   - integer 1/2 - default: 1
%                  Allows to choice the set of coefficients the method
%                  1 - RDE43L a 6-level L-stable method of order 4 and 
%                       4 function evaluations at each step.
%                  2 - RDE43S a 6-level A(x)-stable method (x > 87 deg)
%                       of order 4 and 6 function evaluations at each step.
%
%   In a future version of this code the user will be able to set in OPTIONS
%   many other parameters set as default in this code (see [3 sect.5]).

%   BIBLIOGRAPHY:
%   
%   [1] M. Hochbruck, C. Lubich, and H Selhofer. Exponential integrators for large systems of 
%       differential equations, SIAM J. Sci. Comput. 19 (1998), 1552-1574.
%   [2] I. Moret, On RD-rational Krylov approximations to the core-functions of exponential 
%       integrators, preprint 2005.
%   [3] P. Novati, On the construction of Restricted-Denominator Exponential W-methods, preprint
%       available at http://univaq.it/~novati.
%
%
%   Paolo Novati, 24-Apr-2006.

ct = cputime;
true = 1; false = ~true;

% Initialization

if nargin < 4
    options = [];
    if nargin < 3
        error('not enough input arguments')
    end
end

rtol = odeget(options,'RelTol',1e-3);
atol = odeget(options,'AbsTol',1e-6);
nctr = odeget(options,'NormControl','on');
statson = odeget(options,'Stats','on');
instp = odeget(options,'InitialStep',[]);
hmax = odeget(options,'MaxStep',[]);
if isfield(options,'JPattern')
    Z=options.JPattern;
else
    Z=[];
end
if isfield(options,'Method')
    method = options.Method;
else
    method = 1;
end

if statson
    stats = struct(...
        'nfevals',0,...
        'nsteps',0,...
        'nrejected',0,...
        'npds',0,...
        'ndecomps',0,...
        'nsolves',0,...
        'meankry',0,...
        'CPUtime',0);
end

neq = length(y0);
threshold = atol/rtol;
jthresh = 1e-12;

t0 = tspan(1); tf = tspan(2); 
T = t0; t(1) = t0; y(:,1) = y0; 
htspan = tf - t0;
hmin = 16 * eps * abs(t0);
if isempty(hmax)
    hmax = htspan / 1;
end

if isempty(instp)
    % Compute an initial step using y'(t)
    h = min(hmax, htspan);
    f0 = feval(fode,t0,y0);
    stats.nfevals = stats.nfevals + 1;
    if strcmp(nctr,'on')
        rh = (norm(f0) / max(norm(y0),threshold)) / (0.9 * rtol^(1/2));
    else
        rh = norm(f0 ./ max(abs(y0),threshold),inf) / (0.8 * rtol^(1/3));
    end
    if h * rh > 1
        h = 1 / rh;
    end
    h = max(h, hmin);
else
    h = min(hmax, max(hmin, instp));
end

% Method Parameters

[ga,A,G,b,e,alfa,nfev] = coeff(method);

% Scaling factors for the RD-rational method

sfbase = 5; sfmax = 3; sfmin = 1/3; sf = sfbase;
str = ['(exp(' num2str(sf) '*(1-1./x))-1)./(' num2str(sf) '*(1-1./x))']; fun = inline(str);


I=speye(neq);
if isempty(Z)
    [An,fac,gg,nF] = numjac(fode,t0,y0,f0,jthresh,[],0);
    An = sparse(An); Z = spones(An);
else
    [An,fac,gg,nF] = numjac(fode,t0,y0,f0,jthresh,[],0,Z,[]);
end
[L,U] = lu(I - ga * h / sf * An);
stats.nfevals = stats.nfevals + nF;
stats.npds = stats.npds + 1;
stats.ndecomps = stats.ndecomps + 1; 

% The main loop
f(:,1) = f0;
n = 1; oldh1 = h; oldh2 = h;
done = false;
seqwj = 0; olderr=1;
while ~done
    hmin = 16 * eps * abs(T);
    h = max(hmin, h); 
    if h >= abs(tf - T)
        h = tf - T;
        done = true;
    end

    % Loop for advancing one step  
    nofailed = false;
    while ~nofailed
        w = min(max(abs(y(:,n)), threshold)); 
        
        if done
            nit = 0;
            sf = sfbase;
            [L,U] = lu(I - ga * h / sf * An);
            stats.ndecomps = stats.ndecomps + 1;
            stg(:,1) = f(:,1);
            [lev,nstp] = RDK(L,U,fun,f(:,1),w,sf,h,rtol);
            K(:,1) = lev; nit=nit+nstp;

            for j = 2 : 6
                if j > nfev
                    f(:,j) = f(:,j - 1);
                else
                    f(:,j) = feval(fode,T + alfa(j) * h,y(:,n) + h * (K(:,1:j-1) * A(j,1:j-1)')); 
                end
                stg(:,j) = f(:,j) + h * (An * ((K(:,1:j-1) * G(j,1:j-1)')));
                [lev,nstp] = RDK(L,U,fun,stg(:,j),w,sf,h,rtol);
                K(:,j) = lev; nit = nit + nstp;
            end
            
            stats.nfevals = stats.nfevals + nfev - 1;
            stats.nsolves = stats.nsolves + nit;
            stats.nsteps = stats.nsteps + 1;

            t(n + 1) = tf;
            y(:,n + 1) = y(:,n) + h * K * b';
            break
        end
        
        nit = 0;
        stg(:,1) = f(:,1);
        [lev,nstp] = RDK(L,U,fun,f(:,1),w,sf,h,rtol);
        sfbase = nstp; 
        K(:,1) = lev; nit = nit + nstp;
        for j = 2 : 6
            if j > nfev
                f(:,j) = f(:,j-1);
            else
                f(:,j) = feval(fode,T + alfa(j) * h,y(:,n) + h * (K(:,1:j-1) * A(j,1:j-1)')); 
            end
            stg(:,j) = f(:,j) + h * (An * ((K(:,1:j-1) * G(j,1:j-1)')));
            [lev,nstp] = RDK(L,U,fun,stg(:,j),w,sf,h,rtol);
            K(:,j) = lev; nit = nit + nstp;
        end

        stats.nfevals = stats.nfevals + nfev - 1;
        stats.nsolves = stats.nsolves + nit;

        errv = h * K * e';
        
        % Weighted error estimate
        if strcmp(nctr,'on')
            err = norm( errv ./(rtol * y(:,n) + atol)) * neq^(-1/2);
        else
            err = norm( errv ./ max(abs(y(:,n)),threshold),inf); 
        end

        if err <= 1 % accepted step
            
            seqwj = seqwj + 1;
            nofailed = true;
            T = T + h; t(n + 1) = T; 
            y(:,n + 1) = y(:,n) + h * K * b';
            stats.nsteps = stats.nsteps + 1;
            fac = gus_fact(err,olderr,oldh1,oldh2);
            h = min(oldh1 * min(2, max(0.5 , 0.8 * fac)), hmax); 
            n = n + 1;
            
            f(:,1) = feval(fode,T,y(:,n));
            stats.nfevals = stats.nfevals + 1;
            if abs(h - oldh1) / oldh1 > 0.3 | seqwj > 30
                seqwj = 0;
                [An,fac,gg,nF] = numjac(fode,T,y(:,n),f(:,1),jthresh,[],0,Z,gg);
                sf = sfbase;
                [L,U] = lu(I - ga * h / sf * An);
                stats.nfevals = stats.nfevals + nF;
                stats.npds = stats.npds + 1;
                stats.ndecomps = stats.ndecomps + 1;
            else
                sf = sf * h / oldh1;
                if sf > sfbase * sfmax | sf < sfbase * sfmin
                    sf = sfbase;
                    [L,U] = lu(I - ga * h / sf * An);
                    stats.ndecomps = stats.ndecomps + 1;
                end
            end
            str = ['(exp(' num2str(sf) '*(1-1./x))-1)./(' num2str(sf) '*(1-1./x))']; fun = inline(str);
            oldh2= oldh1;oldh1 = h;  olderr = err;
        else % rejected step
            fac = gus_fact(err,olderr,oldh1,oldh2);
            h = min(oldh1 * min(2, max(0.5 , 0.8 * fac)), hmax);
            
            sf = sf * h / oldh1;
            if sf > sfbase * sfmax | sf < sfbase * sfmin
                sf = sfbase;
                [L,U] = lu(I - ga * h / sf * An);
                stats.ndecomps = stats.ndecomps + 1;
            end
            str = ['(exp(' num2str(sf) '*(1-1./x))-1)./(' num2str(sf) '*(1-1./x))']; fun = inline(str);
            stats.nrejected = stats.nrejected + 1;
        end
    end

end

y = y.';

if statson                   % print cost statistics
    stats.CPUtime = cputime - ct;
    stats.meankry = stats.nsolves / 6 / (stats.nsteps + stats.nrejected);
    fprintf('%g successful steps\n', stats.nsteps);
    fprintf('%g failed attempts\n', stats.nrejected);
    fprintf('%g function evaluations\n', stats.nfevals);
    fprintf('%g partial derivatives\n', stats.npds);
    fprintf('%g LU decomposition\n', stats.ndecomps);
    fprintf('%g solutions of linear systems\n', stats.nsolves);
    fprintf('%g mean number of Krylov iterations\n', stats.meankry);
    fprintf('%g CPU time (seconds)\n', stats.CPUtime);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Krylov iteration %%%%%%%%%%%%%%%%%%%%

function [u,nstp] = RDK(L,U,fun,g,w,sf,h,rtol);

N = length(L);
nr = norm(g);
v(:,1) = g / nr;

PH = 1;
for k = 1 : 10
    
   nstp = k;
   proj = zeros(N,1);
   z = U \ (L \ v(:,k));
   for j = 1 : k
      H(j,k) = v(:,j)' * z;
      proj = proj + H(j,k) * v(:,j);	
   end
   vv = z - proj;
   H(k+1,k) = norm(vv); PH = PH * H(k + 1,k);
   errn = nr * PH * exp(sf - k) * 2^(k + 2) * k^k / sf^(k + 1);
   if h * errn / w < rtol
       break
   end
   v(:,k+1) = vv / H(k + 1,k);
end
tmp1 = funm2(H(1:nstp,:),fun);
u = v(:,1:nstp) * (nr * tmp1(:,1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%  Classocal Factor for the stepsize selection %%%%%%%%%%%%%%

%function fac = cla_fact(err);

%fac = err^(-1/4);

%%%%%%%%%%%%%%%%%%%%%%%%%%% Gustaffson Factor %%%%%%%%%%%%%%%%%%%%%%%%%%%

function fac = gus_fact(err,olderr,h,oldh);

fac = (h/oldh)*(olderr/(err^2))^(1/4);;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  Methods  %%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ga,A,G,b,e,alfa,nfev] = coeff(method);

switch method
    
case 1 % RDE43L
    
    ga=3.333333333333333e-001;
    
    A = zeros(6,6);
    A(2,1) = 1.273220037500351e-001;
    A(3,1) = -9.226441186718157e-001;
    A(3,2) = 1.422644118671816e+000;
    A(4,1) = 3.000000000000000e-001;
    A(5,1) = A(4,1); A(6,1) = A(5,1);
    A(4,2) = 2.000000000000000e-001;
    A(5,2) = A(4,2); A(6,2) = A(5,2);
    A(4,3) = 5.000000000000000e-001;
    A(5,3) = A(4,3); A(6,3) = A(5,3); 
    
    G=zeros(6,6);
    G(2,1) = -9.803381697176562e-002;
    G(3,1) = 8.275955106271749e-001;
    G(3,2) = -6.726441186718158e-001;
    G(4,1) = 8.382963781060364e-001;
    G(5,1) = 1.606553370833683e-001;
    G(6,1) = 7.893003817339626e-001;
    G(4,2) = -2.000000000000000e-001;
    G(5,2) = G(4,2); G(6,2) = G(4,2);
    G(4,3) = -5.000000000000000e-001;
    G(5,3) = G(4,3); G(6,3) = G(5,3);
    G(5,4) = 3.934466291663161e-002;
    
    b(1) = 1.666666666666667e-001;
    b(2) = 0;
    b(3) = 6.666666666666667e-001;
    b(4) = 2.000000000000000e-001;
    b(5) = 5.000000000000000e-001;
    b(6) = -5.333333333333333e-001;
    
    e(1) = -5.802734371614282e-002;
    e(2) = 8.921028305556701e-002;
    e(3) = -3.964901469136306e-002;
    e(4) = 5.223527511893777e-001;
    e(5) = 1.944665749589447e-002;
    e(6) = -5.333333333333333e-001;
    
    alfa(1) = 0;
    alfa(2) = 1.273220037500351e-001;
    alfa(3) = 5.000000000000000e-001;
    alfa(4) = 1;
    alfa(5) = 1;
    alfa(6) = 1;
    
    nfev = 4;
    
case 2 % RDE43S
    
    ga = 2.300000000000000e-001;
    
    A = zeros(6,6);
    A(2,1) = 5.000000000000000e-001;
    A(3,1) = 1.807491994894428e+001;
    A(3,2) = -1.727491994894428e+001;
    A(4,1) = 1.447619738931182e+001;
    A(4,2) = -1.363059573356344e+001;
    A(4,3) = 5.439834425162145e-002;
    A(5,1) = 2.000000000000000e-001; 
    A(5,2) = 5.000000000000000e-001; 
    A(5,3) = 6.000000000000000e-001;  
    A(5,4) = -3.000000000000000e-001;
    A(6,1) = 1.684259259259259259e-001;
    A(6,2) = 6.455555555555555556e-001;
    A(6,3) = 1.319444444444444444e-001;
    A(6,4) = -1.759259259259259259e-001; 
    A(6,5) = 2.300000000000000e-001;
    
    G = zeros(6,6);
    G(2,1) = -5.896681739130403e-001;
    G(3,1) = -2.184086875983272e+001;
    G(3,2) = 1.777491994894428e+001;
    G(4,1) = -1.824368347564747e+001;
    G(4,2) = 1.338336415461607e+001;
    G(4,3) = 7.533995890475487e-005;
    G(5,1) = -3.157407407407407407e-002;
    G(5,2) = 1.455555555555555556e-001; 
    G(5,3) = -4.680555555555555556e-001; 
    G(5,4) = 1.2407407407407407407e-001;
    G(6,5) = -2.300000000000000e-001;
    
    b(1) = 1.684259259259259e-001;
    b(2) = 6.455555555555556e-001;
    b(3) = 1.319444444444445e-001;
    b(4) = -1.759259259259259e-001;
    b(5) = 0;
    b(6) = 2.300000000000000e-001;
    
    e(1) = 0;
    e(2) = 0;
    e(3) = 0;
    e(4) = 0;
    e(5) = -2.300000000000000e-001;
    e(6) = 2.300000000000000e-001;
    
    alfa(1) = 0;
    alfa(2) = .5;
    alfa(3) = .8;
    alfa(4) = .9;
    alfa(5) = 1;
    alfa(6) = 1;
    
    nfev = 6;
    
end