[idegbr,idegcr,m,a,b,c,y,tcos,d,w]=trix(idegbr,idegcr,m,a,b,c,y,tcos,d,w);

function [idegbr,idegcr,m,a,b,c,y,tcos,d,w]=trix(idegbr,idegcr,m,a,b,c,y,tcos,d,w);
persistent i ip k kb kc l lint mm1 x xx z ; 

if isempty(x), x=0; end;
if isempty(xx), xx=0; end;
if isempty(z), z=0; end;
if isempty(i), i=0; end;
if isempty(ip), ip=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(lint), lint=0; end;
if isempty(mm1), mm1=0; end;
%***BEGIN PROLOGUE  TRIX
%***SUBSIDIARY
%***PURPOSE  Subsidiary to GENBUN
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (TRIX-S, CMPTRX-C)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     subroutine to solve a system of linear equations where the
%     coefficient matrix is a rational function in the matrix given by
%     TRIDIAGONAL  ( . . . , A(I), B(I), C(I), . . . ).
%
%***SEE ALSO  GENBUN
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  TRIX
%
a_shape=size(a);a=reshape(a,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
y_shape=size(y);y=reshape(y,1,[]);
tcos_shape=size(tcos);tcos=reshape(tcos,1,[]);
d_shape=size(d);d=reshape(d,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
if isempty(kb), kb=0; end;
if isempty(kc), kc=0; end;
%***FIRST EXECUTABLE STATEMENT  TRIX
mm1 = fix(m - 1);
kb = fix(idegbr + 1);
kc = fix(idegcr + 1);
l =fix(fix((idegbr+1)./(idegcr+1)));
lint = 1;
for k = 1 : idegbr;
x = tcos(k);
if( k==l )
i = fix(idegbr + lint);
xx = x - tcos(i);
for i = 1 : m;
w(i) = y(i);
y(i) = xx.*y(i);
end; i = fix(m+1);
end;
z = 1../(b(1)-x);
d(1) = c(1).*z;
y(1) = y(1).*z;
for i = 2 : mm1;
z = 1../(b(i)-x-a(i).*d(i-1));
d(i) = c(i).*z;
y(i) =(y(i)-a(i).*y(i-1)).*z;
end; i = fix(mm1+1);
z = b(m) - x - a(m).*d(mm1);
if( z~=0. )
y(m) =(y(m)-a(m).*y(mm1))./z;
else;
y(m) = 0.;
end;
for ip = 1 : mm1;
i = fix(m - ip);
y(i) = y(i) - d(i).*y(i+1);
end; ip = fix(mm1+1);
if( k==l )
for i = 1 : m;
y(i) = y(i) + w(i);
end; i = fix(m+1);
lint = fix(lint + 1);
l =fix(fix((lint.*kb)./kc));
end;
end; k = fix(idegbr+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
tcos_shape=zeros(tcos_shape);tcos_shape(:)=tcos(1:numel(tcos_shape));tcos=tcos_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
end
%DECK TSTURM