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

function [wm,iwm,x,tem]=slvs(wm,iwm,x,tem);
%***BEGIN PROLOGUE  SLVS
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DEBDF
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (SLVS-S, DSLVS-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%   SLVS solves the linear system in the iteration scheme for the
%   integrator package DEBDF.
%
%***SEE ALSO  DEBDF
%***ROUTINES CALLED  SGBSL, SGESL
%***COMMON BLOCKS    DEBDF1
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   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  SLVS
%
%LLL. OPTIMIZE
persistent di hl0 i meband ml mu phl0 r ; 

if isempty(i), i=0; end;
global debdf1_12; if isempty(debdf1_12), debdf1_12=0; end;
global debdf1_10; if isempty(debdf1_10), debdf1_10=zeros(1,14); end;
global debdf1_11; if isempty(debdf1_11), debdf1_11=zeros(1,6); end;
global debdf1_13; if isempty(debdf1_13), debdf1_13=0; end;
global debdf1_14; if isempty(debdf1_14), debdf1_14=0; end;
global debdf1_15; if isempty(debdf1_15), debdf1_15=0; end;
global debdf1_18; if isempty(debdf1_18), debdf1_18=0; end;
if isempty(meband), meband=0; end;
global debdf1_16; if isempty(debdf1_16), debdf1_16=0; end;
global debdf1_17; if isempty(debdf1_17), debdf1_17=0; end;
if isempty(ml), ml=0; end;
if isempty(mu), mu=0; end;
global debdf1_19; if isempty(debdf1_19), debdf1_19=0; end;
global debdf1_22; if isempty(debdf1_22), debdf1_22=0; end;
global debdf1_23; if isempty(debdf1_23), debdf1_23=0; end;
global debdf1_20; if isempty(debdf1_20), debdf1_20=0; end;
global debdf1_24; if isempty(debdf1_24), debdf1_24=0; end;
global debdf1_21; if isempty(debdf1_21), debdf1_21=0; end;
global debdf1_1; if isempty(debdf1_1), debdf1_1=0; end;
global debdf1_2; if isempty(debdf1_2), debdf1_2=zeros(1,210); end;
global debdf1_3; if isempty(debdf1_3), debdf1_3=0; end;
global debdf1_4; if isempty(debdf1_4), debdf1_4=0; end;
global debdf1_5; if isempty(debdf1_5), debdf1_5=0; end;
global debdf1_6; if isempty(debdf1_6), debdf1_6=0; end;
global debdf1_7; if isempty(debdf1_7), debdf1_7=0; end;
global debdf1_8; if isempty(debdf1_8), debdf1_8=0; end;
global debdf1_9; if isempty(debdf1_9), debdf1_9=0; end;
if isempty(di), di=0; end;
if isempty(hl0), hl0=0; end;
if isempty(phl0), phl0=0; end;
if isempty(r), r=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 /debdf1/ 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 /debdf1/ debdf1_1 , debdf1_2(210) , debdf1_3 , debdf1_4 , debdf1_5 , debdf1_6 , debdf1_7 ,debdf1_8 , debdf1_9 , debdf1_10(14) , debdf1_11(6) , debdf1_12 ,debdf1_13 , debdf1_14 , debdf1_15 , debdf1_16 , debdf1_17 , debdf1_18 , debdf1_19 ,debdf1_20 , debdf1_21 , debdf1_22 , debdf1_23 , debdf1_24;
%-----------------------------------------------------------------------
% THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM
% A CHORD ITERATION.  IT IS CALLED BY STOD  IF MITER .NE. 0.
% IF MITER IS 1 OR 2, IT CALLS SGESL 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 SGBSL.
% COMMUNICATION WITH SLVS USES THE FOLLOWING VARIABLES..
% WM  = REAL 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.
%-----------------------------------------------------------------------
%***FIRST EXECUTABLE STATEMENT  SLVS
debdf1_12 = 0;
if( debdf1_17==3 )
%
phl0 = wm(2);
hl0 = debdf1_4.*debdf1_3;
wm(2) = hl0;
if( hl0~=phl0 )
r = hl0./phl0;
for i = 1 : debdf1_19;
di = 1.0e0 - r.*(1.0e0-1.0e0./wm(i+2));
if( abs(di)==0.0e0 )
debdf1_12 = -1;
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.0e0./di;
end;
end; i = fix(debdf1_19+1);
end;
for i = 1 : debdf1_19;
x(i) = wm(i+2).*x(i);
end; i = fix(debdf1_19+1);
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;
elseif( debdf1_17==4 || debdf1_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,debdf1_19,ml,mu,iwm(sub2ind(size(iwm),max(21,1)):end),x]=sgbsl(wm(sub2ind(size(wm),max(3,1)):end),meband,debdf1_19,ml,mu,iwm(sub2ind(size(iwm),max(21,1)):end),x,0);
else;
debdf1_19_orig=debdf1_19;    [wm(sub2ind(size(wm),max(3,1)):end),debdf1_19,dumvar3,iwm(sub2ind(size(iwm),max(21,1)):end),x]=sgesl(wm(sub2ind(size(wm),max(3,1)):end),debdf1_19,debdf1_19,iwm(sub2ind(size(iwm),max(21,1)):end),x,0);    debdf1_19(dumvar3~=debdf1_19_orig)=dumvar3(dumvar3~=debdf1_19_orig);
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;
end;
%----------------------- end OF SUBROUTINE SLVS -----------------------
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 SMOUT