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

%   Reserch code
%
%   ROSWB solves stiff ordinary differential equations by ROW-methods
%
%   [t,y] = ROSWB(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. Even if the syntax is ODEFUN(T,Y), 
%   ODEFUN must be independent of T.
%
%   [t,y] = ROSWB(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
%   integration properties replaced by values in OPTIONS, an argument created
%   with the ODESET function. See ODESET for details. ROSWB supports the
%   most important features of the matlab ODE solver. Moreover in this code
%   there are a few other options available:
%
%       'Method'   - integer between 1 and 10 - default: 1
%                  Allows to choice the set of coefficients of the ROW-method used
%                  1 - WB34 [2] a 6-level stiffly-accurate method of order 4
%                  2 - GRK4T [3] a 4-level strongly A-stable method of order 4
%                  3 - GRK4A [3] a 4-level strongly A-stable method of order 4
%                  4 - VELDD [4] a 4-level strongly A-stable method of order 4
%                  5 - LSTAB [1] a 4-level L-stable method of order 4
%                  6 - RODAS [1] a 6-level stiffly-accurate method of order 4
%                  7 - W23B [2] a 4-level stiffly-accurate method of order 3
%                  8 - ROS3P [5] a 3-level method of order 3
%                  9 - ROS3PW [6] a 3-level method of order 3
%                  10 - RODASP [5] a 6-level method of order 4
%       'BroydenUpdate' - integer between 0 and 4 - default: 0
%                  Enables the choice of the rank-1 update of the Jacobian
%                  0 - no update, the Jacobian is computed at each step
%                  1 - Broyden's good update
%                  2 - Broyden's bad update
%                  3 - Schubert's update (sparse update)
%                  4 - Constant Jacobian 
%       'Recover'  - 'on'/'off' - default: 'off'
%                  Allows to recover the Jacobian when a failed attempt occures
%                  using the 'BroydenUpdate' 1,2,3,4. 
%
%   Since the use of these additional options is not allowed by ODESET, OPTIONS must be
%   createted using 'struct'. Example: 
%
%   options = struct('RelTol',1e-9,'Method',10,'BroydenUpdate',0,'Recover','off')
%
%   solve the ODE with RODASP, with RelTol = 1e-9 and AbsTol = 1e-6 (by default)
%
%   REMARK: For large problems, if the Jacobian is not available, it is
%   strongly recommended to give at least its sparsity pattern. Actually,
%   if not given, after the first call to numjac the code automatically generates it, but
%   in some cases the methods does not work well.
%
%   BIBLIOGRAPHY:
%   
%   [1] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II, second ed., 
%       Springer, Berlin, 1996.
%   [2] P. Novati, Some secant approximations for Rosenbrock W-methods,
%       preprint available at http://univaq.it/~novati.
%   [3] P. Kaps, P. Rentrop, Generalized RungeKutta methods of order four with stepsize control 
%       for stiff ordinary differential equations, Numer. Math. 38 (1979) 55-68.
%   [4] M.V. van Veldhuizen, D-stability and Kaps-Rentrop methods, Computing 32 (1984), 229-237.
%   [5] J. Lang, J.G. Verwer, ROS3P - An accurate third-order Rosenbrock solver
%       designed for parabolic problems, BIT 41 (2001), 731-738.
%   [6] J. Rang, A. Angermann, New Rosenbrock W-methods of order 3 for PDAEs
%       of index 1, BIT 43 (2003), 1-23.
%   [7] G. Steinebach, Order-reduction of ROW-methods for DAEs and method
%       of lines applications, Preprint-Nr. 1741, FB Mathematik, TH Darmstadt (1995).
%
%
%   Paolo Novati, 21-Nov-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);
statson = odeget(options,'Stats','off');
instp = odeget(options,'InitialStep',[]);
hmax = odeget(options,'MaxStep',[]);
fjac = odeget(options,'Jacobian',[]);
Z = odeget(options,'JPattern',[]);
if statson
    stats = struct(...
        'nfevals',0,...
        'nsteps',0,...
        'nrejected',0,...
        'npds',0,...
        'ndecomps',0,...
        'nsolves',0,...
        'CPUtime',0);
end

if isfield(options,'Method')
    method = options.Method;
else
    method = 1;
end
if isfield(options,'BroydenUpdate')
    broyden = options.BroydenUpdate;
else
    broyden = 0;
end
if isfield(options,'Recover')
    if broyden > 0
        recover = options.Recover;
    else
        recover = 'off';
    end
else
    recover = 'off';
end

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

% Time settings
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/10;
end

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

f(:,1) = f0;
I = speye(neq);

% Handle the Jacobian
if isempty(fjac)
    if isempty(Z)
        [An,fc,gg,nF] = numjac(fode,t0,y0,f0,jthresh,[],0);
        An = sparse(An); Z = spones(An);
        stats.nfevals = stats.nfevals + nF;
    else
        [An,fc,gg,nF] = numjac(fode,t0,y0,f0,jthresh,[],0,Z,[]);
        An=sparse(An);
        stats.nfevals = stats.nfevals + nF;
    end
else
    An=feval(fjac,t0,y0);
end
stats.npds = stats.npds + 1;

% Method Parameters
[ga,a,c,m,e,alfa,nfev,nlev,ord] = coeff(method);

if broyden > 0
    facmax = 2; facmin = 0.2; facs = 0.75; ord = ord - 1;
else
    facmax = 5; facmin = 0.2; facs = 0.75;
end

[L,U] = lu(I - h * ga * An);
if broyden == 1 | broyden == 2
    vx = 0; vy = 0; avx = 0; avy = 0; sk = [];
end
stats.ndecomps = stats.ndecomps + 1;

% The main loop
n = 1;
done = false; fouls = 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
        
        if broyden == 0 | broyden >=3 
            u(:,1) = h * ga * (U \ (L \ (f(:,1)))); 
            for j = 2 : nlev
                if j > nfev
                    f(:,j) = f(:,j-1);
                else
                    f(:,j) = feval(fode, T + alfa(j) * h, y(:,n) + (u(:,1:j-1) * a(j,1:j-1)')); 
                end
                u(:,j) = h * ga * (U \ (L \ ((f(:,j) + 1/h * (u(:,1:j-1) * c(j,1:j-1)')))));
            end
            
        else
            
            u(:,1) = h * ga * (U \ (L \ f(:,1)) + vx * (vy' * f(:,1))); 
            for j = 2 : nlev
                if j > nfev
                    f(:,j) = f (:,j-1);
                else
                    f(:,j) = feval(fode, T + alfa(j) * h, y(:,n) + (u(:,1:j-1) * a(j,1:j-1)')); 
                end
                tmp = f(:,j) + 1/h * u(:,1:j-1) * c(j,1:j-1)';
                u(:,j) = h * ga * (U \ (L \ tmp) + vx * (vy' * tmp));
            end
        end
        stats.nfevals = stats.nfevals + nfev - 1;
        stats.nsolves = stats.nsolves + nlev;

        if done
            T = T + h; t(n + 1) = T; 
            y(:,n + 1) = y(:,n) + u * m';
            stats.nsteps = stats.nsteps + 1;
            break
        end            
                
        errv = u * e';
        
        % Weighted error estimate
        err = norm( errv ./(rtol * y(:,n) + atol)) * neq^(-1/2) * rtol;
        if err <= rtol 
            fouls = 0;
            nofailed = true;
            T = T + h; t(n + 1) = T; 
            y(:,n + 1) = y(:,n) + u * m';
            stats.nsteps = stats.nsteps + 1;
            fac = cla_fact(rtol,err,h,ord);
            h = min(h * min(facmax, max(facmin , facs * fac)), hmax); 
            n = n + 1;
            fo = f(:,1);
            f(:,1) = feval(fode,T,y(:,n));
            stats.nfevals = stats.nfevals + 1;
            
            if broyden < 4
                if broyden > 0
                    dy = y(:,n) - y(:,n-1); 
                    df = f(:,1) - fo;
                    switch broyden
                    case 1 % Good Broyden's update
                        [avx,avy,sk,vx,vy] = update_jacobian_g(An,L,U,dy,df,h,oldh,ga,sk,avx,avy,vx,vy);
                        stats.nsolves = stats.nsolves + 2;
                    case 2 % Bad Broyden's update
                        [vx,vy] = update_jacobian_b(L,U,dy,df,h,ga,vx,vy);
                        stats.nsolves = stats.nsolves + 1;
                    case 3 % Schubert's update
                        An = update_jacobian_s(An,dy,df,Z,neq);
                        [L,U] = lu(I - h * ga * An);
                        stats.ndecomps = stats.ndecomps + 1;
                    end
                else
                    if isempty(fjac)
                        [An,fc,gg,nF] = numjac(fode,T,y(:,n),f(:,1),jthresh,fc,0,Z,gg);
                        stats.nfevals = stats.nfevals + nF;
                    else
                        An = feval(fjac,T,y(:,n));
                    end
                    stats.npds = stats.npds + 1;
                    [L,U] = lu(I - h * ga * An);
                    stats.ndecomps = stats.ndecomps + 1;
                end
            else
                [L,U] = lu(I - h * ga * An);
                stats.ndecomps = stats.ndecomps + 1;
            end
            oldh = h;
        else % failed attempts
            fac = cla_fact(rtol,err,h,ord);
            h = min(h * min(facmax, max(facmin , facs * fac)), hmax);
            stats.nrejected = stats.nrejected + 1;
            if strcmp(recover,'on') & broyden > 0 & fouls == 0
                if isempty(fjac)
                    [An,fc,gg,nF] = numjac(fode,T,y(:,n),f(:,1),jthresh,[],0,Z,gg);
                    
                    stats.nfevals = stats.nfevals + nF;
                else
                    An = feval(fjac,T,y(:,n));
                end
                stats.npds = stats.npds + 1;
            end
            [L,U] = lu(I - h * ga * An);
            if broyden == 1 | broyden ==2
                vx = 0; vy = 0; avx = 0; avy = 0; sk = [];
            end
            stats.ndecomps = stats.ndecomps + 1;
            fouls = fouls + 1;
        end
    end

end

y = y.';

if strcmp(statson,'on')                   % print cost statistics
    stats.CPUtime = cputime - ct;
    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 CPU time (seconds)\n', stats.CPUtime);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Good Broyden's update %%%%%%%%%%%%%%%%%%%%%

function [navx,navy,sk,nvx,nvy] = update_jacobian_g(A,L,U,dy,df,h,oldh,ga,sk,avx,avy,vx,vy);

nn = length(sk);

if nn == 0
    Ady = A*dy;
else
    Ady = prod(sk) * (A * dy);
    for k = 1 : nn
        Ady = Ady + prod(sk(k:nn)) * (avx(:,k) * (avy(:,k)' * dy));
    end
end
vk = (df * h - Ady * oldh) / (oldh * (dy' * dy));
sk = [sk oldh/h];

if vx == 0
    navx = vk; navy = dy;
else
    navx = [avx vk]; navy = [avy dy];
end

tmp = (ga * oldh) * vk;
ck = - (U \ (L\ tmp) + vx * (vy' * tmp));

if vx == 0
    nvx = -(ck / (1 + dy' * ck)); nvy = L' \ (U' \ dy) + vy * (vx' * dy);
else
    nvx = [vx - (ck / (1 + dy' * ck))]; nvy = [vy L' \ (U' \ dy) + vy * (vx' * dy)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bad Broyden's update %%%%%%%%%%%%%%%%%%%%%%

function [nvx,nvy] = update_jacobian_b(L,U,dy,df,h,ga,vx,vy);


vk = dy - h * ga * df;
vn = vk' / (vk' * vk);
if vx == 0
    nvy = vn'; nvx = dy - (U \ (L \ vk));
else
    nvx = [vx dy - (U \ (L \ vk)) - vx * (vy' * vk)]; nvy = [vy vn'];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Schubert's update %%%%%%%%%%%%%%%%%%%%%%%%%

function An = update_jacobian_s(A,dy,df,Z,n);

S = spdiags(dy,0,n,n) * Z'; 
k = find(Z); [i,j] = ind2sub([n n],k);
vect = diag(S' * S); vect(find(abs(vect) < 10e-9)) = 1;
tmp = (df - A * dy) ./ vect;
Z(k) = tmp(i) .* dy(j);
An = A + Z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  Classical factor %%%%%%%%%%%%%%%%%%%%%%%%

function fac = cla_fact(rtol,err,h,ord);

fac = (rtol/err)^(1/ord);

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

function [ga,a,c,m,e,alfa,nfev,nlev,ord] = coeff(method);

switch method
case 1 % W34B
    
    ga = 5.728160624821350e-001;
    
    a(2,1) = 9.077957726023402e-001;
    a(3,1) = 1.488057317686709e+000;
    a(4,1) = 5.706683153980792e-001;
    a(5,1) = 1.777502264129759e+000;
    a(6,1) = 1.777502264129759e+000;
    a(3,2) = 1.090896667173186e+000;
    a(4,2) = -2.742670573870741e+000;
    a(5,2) = 2.760496374413255e+000;
    a(6,2) = 2.760496374413255e+000;
    a(4,3) = 6.826575262514714e-001;
    a(5,3) = -5.548392770606871e-001;
    a(6,3) = -5.548392770606871e-001;
    a(5,4) = -2.755429936381875e-001;
    a(6,4) = -2.755429936381875e-001;
    a(6,5) = 1;
    
    c(2,1) = -1.584794547605153e+000;
    c(3,1) = -1.354926096553657e+000;
    c(4,1) = -5.471574707103881e-001;
    c(5,1) = -9.092693797373024e-001;
    c(6,1) = -9.092693797373024e-001;
    c(3,2) = 1.981427137267801e+000;
    c(4,2) = 4.043496562125793e+000;
    c(5,2) = 3.733616763910639e+000;
    c(6,2) = 3.733616763910639e+000;
    c(4,3) = -2.777778234947648e+000;
    c(5,3) = -2.816848938364382e+000;
    c(6,3) = -2.816848938364382e+000;
    c(5,4) = -6.765567623376954e-001;
    c(6,4) = -6.765567623376954e-001;
    c(6,5) = -1.745761101158346e+000;

    m(1) = 1.777502264129759e+000; 
    m(2) = 2.760496374413255e+000;
    m(3) = -5.548392770606871e-001;
    m(4) = -2.755429936381874e-001;
    m(5) = 1;
    m(6) = 1;
    
    e(6) = 1;
    
    alfa(2) = 0.52;
    alfa(3) = 0.91;
    alfa(4) = 3.104974864299539e-001;
    alfa(5) = 1;
    alfa(6) = 1;
    
    nfev = 6;
    nlev = 6;
    ord = 4;
        
case 2 % GRK4T
    
    a(2,1) = 0.2000000000000000e+001;
    a(3,1) = 0.4524708207373116e+001;
    a(3,2) = 0.4163528788597648e+001;
    c(2,1) = -0.5071675338776316e+001;
    c(3,1) = 0.6020152728650786e+001;
    c(3,2) = 0.1597506846727117e+000;
    c(4,1) = -0.1856343618686113e+001;
    c(4,2) = -0.8505380858179826e+001;
    c(4,3) = -0.2084075136023187e+001;
    m(1) = 0.3957503746640777e+001;
    m(2) = 0.4624892388363313e+001;
    m(3) = 0.6174772638750108e+000;
    m(4) = 0.1282612945269037e+001;
    e(1) = 0.2302155402932996e+001;
    e(2) = 0.3073634485392623e+001;
    e(3) = -0.8732808018045032e+000;
    e(4) = -0.1282612945269037e+001;
    ga = 0.2310000000000000e+000;
    alfa(2) = 0.4620000000000000e+000;
    alfa(3) = 0.8802083333333334e+000;
    alfa(4) = 0.8802083333333334e+000;
    
    nfev = 3;
    nlev = 4;
    ord = 4;

case 3 % GRK4A
    
    a(2,1) = 0.1108860759493671e+001;
    a(3,1) = 0.2377085261983360e+001;
    a(3,2) = 0.1850114988899692e+000;
    c(2,1) = -0.4920188402397641e+001;
    c(3,1) = 0.1055588686048583e+001;
    c(3,2) = 0.3351817267668938e+001;
    c(4,1) = 0.3846869007049313e+001;
    c(4,2) = 0.3427109241268180e+001;
    c(4,3) = -0.2162408848753263e+001;
    m(1) = 0.1845683240405840e+001;
    m(2) = 0.1369796894360503e+000;
    m(3) = 0.7129097783291559e+000;
    m(4) = 0.6329113924050632e+000;
    e(1) = 0.4831870177201765e-001;
    e(2) = -0.6471108651049505e+000;
    e(3) = 0.2186876660500240e+000;
    e(4) = -0.6329113924050632e+000;
    ga = 0.3950000000000000e+000;
    alfa(2) = 0.4380000000000000e+000;
    alfa(3) = 0.8700000000000000e+000;
    alfa(4) = 0.8700000000000000e+000;
    
    nfev = 3;
    nlev = 4;
    ord = 4;
    
case 4 % VELDD
    
    a(2,1) = 0.2000000000000000e+01;
    a(3,1) = 0.4812234362695436e+01;
    a(3,2) = 0.4578146956747842e+01;
    c(2,1) = -0.5333333333333331e+01;
    c(3,1) = 0.6100529678848254e+01;
    c(3,2) = 0.1804736797378427e+01;
    c(4,1) = -0.2540515456634749e+01;
    c(4,2) = -0.9443746328915205e+01;
    c(4,3) = -0.1988471753215993e+01;
    m(1) = 0.4289339254654537e+01;
    m(2) = 0.5036098482851414e+01;
    m(3) = 0.6085736420673917e+00;
    m(4) = 0.1355958941201148e+01;
    e(1) = 0.2175672787531755e+01;
    e(2) = 0.2950911222575741e+01;
    e(3) = -0.7859744544887430e+00;
    e(4) = -0.1355958941201148e+01;
    ga = 0.2257081148225682e+00;
    alfa(2) = 0.4514162296451364e+00;
    alfa(3) = 0.8755928946018455e+00;
    alfa(4) = 0.8755928946018455e+00;

    nfev = 3;
    nlev = 4;
    ord = 4;
    
case 5 % LSTAB
    
    a(2,1) = 0.2000000000000000e+01;
    a(3,1) = 0.1867943637803922e+01;
    a(3,2) = 0.2344449711399156e+00;
    c(2,1) = -0.7137615036412310e+01;
    c(3,1) = 0.2580708087951457e+01;
    c(3,2) = 0.6515950076447975e+00;
    c(4,1) = -0.2137148994382534e+01;
    c(4,2) = -0.3214669691237626e+00;
    c(4,3) = -0.6949742501781779e+00;
    m(1) = 0.2255570073418735e+01;
    m(2) = 0.2870493262186792e+00;
    m(3) = 0.4353179431840180e+00;
    m(4) = 0.1093502252409163e+01;
    e(1) = -0.2815431932141155e+00;
    e(2) = -0.7276199124938920e-01;
    e(3) = -0.1082196201495311e+00;
    e(4) = -0.1093502252409163e+01;
    ga = 0.5728200000000000e+00;
    alfa(2) = 0.1145640000000000e+01;
    alfa(3) = 0.6552168638155900e+00;
    alfa(4) = 0.6552168638155900e+00;

    nfev = 3;
    nlev = 4;
    ord = 4;
    
case 6 % ROCOE
    
    alfa(2) = 0.386;
    alfa(3) = 0.21; 
    alfa(4) = 0.63;
    alfa(5) = 1;
    alfa(6) = 1;
        
    m(1) = 1.221224509226641e+000;
    m(2) = 6.019134481288628e+000;
    m(3) = 1.253708332932088e+001;
    m(4) = -6.878860361058932e-001;
    m(5) = 1;
    m(6) = 1;

    e(6) = 1;
    a(2,1) = 0.1544000000000000e+01;
    a(3,1) = 0.9466785280815826e+00;
    a(3,2) = 0.2557011698983284e+00;
    a(4,1) = 0.3314825187068521e+01;
    a(4,2) = 0.2896124015972201e+01;
    a(4,3) = 0.9986419139977817e+00;
    a(5,1) = 0.1221224509226641e+01;
    a(5,2) = 0.6019134481288629e+01;
    a(5,3) = 0.1253708332932087e+02;
    a(5,4) = -0.6878860361058950e+00;
    a(6,1) = a(5,1); 
    a(6,2) = a(5,2);
    a(6,3) = a(5,3);
    a(6,4) = a(5,4);
    a(6,5) = 1;
    c(2,1) = -0.5668800000000000e+01;
    c(3,1) = -0.2430093356833875e+01;
    c(3,2) = -0.2063599157091915e+00;
    c(4,1) = -0.1073529058151375e+00;
    c(4,2) = -0.9594562251023355e+01;
    c(4,3) = -0.2047028614809616e+02;
    c(5,1) = 0.7496443313967647e+01;
    c(5,2) = -0.1024680431464352e+02;
    c(5,3) = -0.3399990352819905e+02;
    c(5,4) = 0.1170890893206160e+02;
    c(6,1) = 0.8083246795921522e+01;
    c(6,2) = -0.7981132988064893e+01;
    c(6,3) = -0.3152159432874371e+02;
    c(6,4) = 0.1631930543123136e+02;
    c(6,5) = -0.6058818238834054e+01;
    ga = 0.25;

    nfev = 6;
    nlev = 6;
    ord = 4;
    
case 7 % W23B
    
    a(2,1) = 1.147140180139521e+000;
    a(3,1) = 2.530586938130458e+000;
    a(3,2) = 1.605996252195329e+000;
    a(4,1) = 2.530586938130458e+000;
    a(4,2) = 1.605996252195329e+000;
    m(1) = 2.530586938130458e+000;
    m(2) = 1.605996252195329e+000;
    m(3) = 0;
    m(4) = 1;
    e(1) = 1.683237534765762e-001;
    e(2) = 9.553251206785370e-001;
    e(3) = -1.089666303771148e-001;
    e(4) = 0.5;
    ga = 4.358665215084590e-001;
    c(2,1) = -2.631861185781065e+000;
    c(3,1) = -1.457265169704310e+000;
    c(3,2) = 1.716687547938875e+000;
    c(4,1) = -1.803139923801317e+000;
    c(4,2) = 8.848030478915137e-001;
    c(4,3) = -1.416993298352020e+000;
    alfa(2) = .5;
    alfa(3) = 1;
    alfa(4) = 1;
    nfev = 3;
    nlev = 4;
    ord = 3;
    
case 8 % ROS3P

    ga = .7886751345948129e+00;
    a(2,1) = 1.267949192431123e+00;
    c(2,1) = -1.607695154586736e+00;
    a(3,1) = 1.267949192431123e+00;
    c(3,1) = -3.464101615137755e+00;
    a(3,2) = 0; 
    c(3,2) = -1.732050807568877e+00;
    alfa(1) = 0; alfa(2) = 1; alfa(3) = 1;
    m(1) = 2;
    e(1)= m(1) - 2.113248654051871e+00;
    m(2) = 5.773502691896258e-01;
    e(2) = m(2) -1;
    m(3) = 4.226497308103742e-01;
    e(3) = m(3) - 4.226497308103742e-01;
    nlev = 3;
    nfev = 2;
    ord = 3;
    
case 9 % ROS3PW
    
    ga = .7886751345948129e+00;
    a(2,1) = 2;
    c(2,1) = -2.53589838486225;
    a(3,1) = 0.63397459621556;
    c(3,1) = -1.62740473580836;
    a(3,2) = 0; 
    c(3,2) = -0.27451905283833;
    alfa(1) = 0; alfa(2) = 2 * ga; alfa(3) = .5;
    m = [1.63397459621556   0.29422863405995   1.07179676972449];
    e = [-0.36047190431930  -0.36047190431930  0];
    nlev = 3;
    nfev = 3;
    ord = 3;
    
case 10 % RODASP
    ga = 0.25;
    a = [0 0 0 0 0 0;
   3.00000000000000 0 0 0 0 0;
   1.83103679348678   0.49551839674339 0 0 0 0;
   2.30437658269264  -0.05249275245744  -1.17679876183277 0 0 0;
  -7.17045496242342  -4.74163667148198 -16.31002631330972  -1.06200404411140 0 0;
  -7.17045496242341  -4.74163667148198 -16.31002631330972  -1.06200404411140 1 0];
        
   c= [ 0 0 0 0 0 0;
 -12  0 0 0 0 0;
  -8.79179517394713  -2.20786558697357  0 0 0 0;
  10.81793056857201   6.78027061142850  19.53485944642405 0 0 0;
  34.19095006749803  15.49671153726025  54.74760875964146  14.16005392148532 0 0;
  34.62605830930666  15.30084976114536  56.99955578662686  18.40807009793090  -5.71428571428571 0];
 
  m = [-7.17045496242341  -4.74163667148197 -16.31002631330972  -1.06200404411140 1 1];
  alfa = [0 0.75 0.21 0.63 1 1];

  e = [0 0 0 0 0 1];
  nlev = 6; nfev = 6; ord = 4;
  
end
