[wm,iwm,x,tem]=dslvs(wm,iwm,x,tem);

function [wm,iwm,x,tem]=dslvs(wm,iwm,x,tem);
%***BEGIN PROLOGUE  DSLVS
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DDEBDF
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (SLVS-S, DSLVS-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%   DSLVS solves the linear system in the iteration scheme for the
%   integrator package DDEBDF.
%
%***SEE ALSO  DDEBDF
%***ROUTINES CALLED  DGBSL, DGESL
%***COMMON BLOCKS    DDEBD1
%***REVISION HISTORY  (YYMMDD)
%   820301  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%   910722  Updated AUTHOR section.  (ALS)
%   920422  Changed DIMENSION statement.  (WRB)
%***end PROLOGUE  DSLVS
%
persistent di hl0 i meband ml mu phl0 r ; 

if isempty(i), i=0; end;
global ddebd1_12; if isempty(ddebd1_12), ddebd1_12=0; end;
global ddebd1_10; if isempty(ddebd1_10), ddebd1_10=zeros(1,14); end;
global ddebd1_11; if isempty(ddebd1_11), ddebd1_11=zeros(1,6); end;
global ddebd1_13; if isempty(ddebd1_13), ddebd1_13=0; end;
global ddebd1_14; if isempty(ddebd1_14), ddebd1_14=0; end;
global ddebd1_15; if isempty(ddebd1_15), ddebd1_15=0; end;
global ddebd1_18; if isempty(ddebd1_18), ddebd1_18=0; end;
if isempty(meband), meband=0; end;
global ddebd1_16; if isempty(ddebd1_16), ddebd1_16=0; end;
global ddebd1_17; if isempty(ddebd1_17), ddebd1_17=0; end;
if isempty(ml), ml=0; end;
if isempty(mu), mu=0; end;
global ddebd1_19; if isempty(ddebd1_19), ddebd1_19=0; end;
global ddebd1_22; if isempty(ddebd1_22), ddebd1_22=0; end;
global ddebd1_23; if isempty(ddebd1_23), ddebd1_23=0; end;
global ddebd1_20; if isempty(ddebd1_20), ddebd1_20=0; end;
global ddebd1_24; if isempty(ddebd1_24), ddebd1_24=0; end;
global ddebd1_21; if isempty(ddebd1_21), ddebd1_21=0; end;
if isempty(di), di=0; end;
global ddebd1_3; if isempty(ddebd1_3), ddebd1_3=0; end;
global ddebd1_4; if isempty(ddebd1_4), ddebd1_4=0; end;
if isempty(hl0), hl0=0; end;
global ddebd1_5; if isempty(ddebd1_5), ddebd1_5=0; end;
global ddebd1_6; if isempty(ddebd1_6), ddebd1_6=0; end;
global ddebd1_7; if isempty(ddebd1_7), ddebd1_7=0; end;
if isempty(phl0), phl0=0; end;
if isempty(r), r=0; end;
global ddebd1_1; if isempty(ddebd1_1), ddebd1_1=0; end;
global ddebd1_2; if isempty(ddebd1_2), ddebd1_2=zeros(1,210); end;
global ddebd1_8; if isempty(ddebd1_8), ddebd1_8=0; end;
global ddebd1_9; if isempty(ddebd1_9), ddebd1_9=0; end;
wm_shape=size(wm);wm=reshape(wm,1,[]);
iwm_shape=size(iwm);iwm=reshape(iwm,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
tem_shape=size(tem);tem=reshape(tem,1,[]);
% common :: ;
%% common /ddebd1/ rownd , rowns(210) , el0 , h , hmin , hmxi , hu ,tn , uround , iownd(14) , iowns(6) , ier ,jstart , kflag , l , meth , miter , maxord , n ,nq , nst , nfe , nje , nqu;
%% common /ddebd1/ ddebd1_1 , ddebd1_2(210) , ddebd1_3 , ddebd1_4 , ddebd1_5 , ddebd1_6 , ddebd1_7 ,ddebd1_8 , ddebd1_9 , ddebd1_10(14) , ddebd1_11(6) , ddebd1_12 ,ddebd1_13 , ddebd1_14 , ddebd1_15 , ddebd1_16 , ddebd1_17 , ddebd1_18 , ddebd1_19 ,ddebd1_20 , ddebd1_21 , ddebd1_22 , ddebd1_23 , ddebd1_24;
%     ------------------------------------------------------------------
%      THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING
%      FROM A CHORD ITERATION.  IT IS CALLED BY DSTOD  IF MITER .NE. 0.
%      IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS.
%      IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
%      MATRIX, AND THEN COMPUTES THE SOLUTION.
%      IF MITER IS 4 OR 5, IT CALLS DGBSL.
%      COMMUNICATION WITH DSLVS USES THE FOLLOWING VARIABLES..
%      WM  = doubleprecision WORK SPACE CONTAINING THE INVERSE DIAGONAL
%      MATRIX IF MITER
%            IS 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE.
%            STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
%            WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
%            WM(1) = SQRT(UROUND) (NOT USED HERE),
%            WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER =
%            3.
%      IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING
%            AT IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS
%            THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS
%            4 OR 5.
%      X   = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION
%            VECTOR ON OUTPUT, OF LENGTH N.
%      TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION.
%      IER = OUTPUT FLAG (IN COMMON).  IER = 0 IF NO TROUBLE OCCURRED.
%            IER = -1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
%      THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
%-----------------------------------------------------------------------
%     BEGIN BLOCK PERMITTING ...EXITS TO 80
%        BEGIN BLOCK PERMITTING ...EXITS TO 60
%***FIRST EXECUTABLE STATEMENT  DSLVS
ddebd1_12 = 0;
if( ddebd1_17==3 )
%
phl0 = wm(2);
hl0 = ddebd1_4.*ddebd1_3;
wm(2) = hl0;
if( hl0~=phl0 )
r = hl0./phl0;
for i = 1 : ddebd1_19;
di = 1.0d0 - r.*(1.0d0-1.0d0./wm(i+2));
%        .........EXIT
if( abs(di)==0.0d0 )
ddebd1_12 = -1;
%     ...EXIT
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
tem_shape=zeros(tem_shape);tem_shape(:)=tem(1:numel(tem_shape));tem=tem_shape;
return;
else;
wm(i+2) = 1.0d0./di;
end;
end; i = fix(ddebd1_19+1);
end;
for i = 1 : ddebd1_19;
x(i) = wm(i+2).*x(i);
%     ......EXIT
end; i = fix(ddebd1_19+1);
elseif( ddebd1_17==4 || ddebd1_17==5 ) ;
%
ml = fix(iwm(1));
mu = fix(iwm(2));
meband = fix(2.*ml + mu + 1);
[wm(sub2ind(size(wm),max(3,1)):end),meband,ddebd1_19,ml,mu,iwm(sub2ind(size(iwm),max(21,1)):end),x]=dgbsl(wm(sub2ind(size(wm),max(3,1)):end),meband,ddebd1_19,ml,mu,iwm(sub2ind(size(iwm),max(21,1)):end),x,0);
else;
%     ......EXIT
ddebd1_19_orig=ddebd1_19;    [wm(sub2ind(size(wm),max(3,1)):end),ddebd1_19,dumvar3,iwm(sub2ind(size(iwm),max(21,1)):end),x]=dgesl(wm(sub2ind(size(wm),max(3,1)):end),ddebd1_19,ddebd1_19,iwm(sub2ind(size(iwm),max(21,1)):end),x,0);    ddebd1_19(dumvar3~=ddebd1_19_orig)=dumvar3(dumvar3~=ddebd1_19_orig);
end;
%     ----------------------- end OF SUBROUTINE DSLVS
%     -----------------------
wm_shape=zeros(wm_shape);wm_shape(:)=wm(1:numel(wm_shape));wm=wm_shape;
iwm_shape=zeros(iwm_shape);iwm_shape(:)=iwm(1:numel(iwm_shape));iwm=iwm_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
tem_shape=zeros(tem_shape);tem_shape(:)=tem(1:numel(tem_shape));tem=tem_shape;
end
%DECK DSMMI2