[m,a,b,c,k,y1,y2,y3,tcos,d,w1,w2,w3]=tri3(m,a,b,c,k,y1,y2,y3,tcos,d,w1,w2,w3);

function [m,a,b,c,k,y1,y2,y3,tcos,d,w1,w2,w3]=tri3(m,a,b,c,k,y1,y2,y3,tcos,d,w1,w2,w3);
persistent i ip k1 k1p1 k2 k2k3k4 k2p1 k3 k3p1 k4 k4p1 kint1 kint2 kint3 l1 l2 l3 lint1 lint2 lint3 mm1 n 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(k1), k1=0; end;
if isempty(k2), k2=0; end;
if isempty(k2k3k4), k2k3k4=0; end;
if isempty(k3), k3=0; end;
if isempty(k4), k4=0; end;
if isempty(kint1), kint1=0; end;
if isempty(kint2), kint2=0; end;
if isempty(kint3), kint3=0; end;
if isempty(l1), l1=0; end;
if isempty(l2), l2=0; end;
if isempty(l3), l3=0; end;
if isempty(lint1), lint1=0; end;
if isempty(lint2), lint2=0; end;
if isempty(lint3), lint3=0; end;
if isempty(mm1), mm1=0; end;
if isempty(n), n=0; end;
%***BEGIN PROLOGUE  TRI3
%***SUBSIDIARY
%***PURPOSE  Subsidiary to GENBUN
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (TRI3-S, CMPTR3-C)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%     subroutine to solve three linear systems whose common 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
%   890206  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  TRI3
a_shape=size(a);a=reshape(a,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
tcos_shape=size(tcos);tcos=reshape(tcos,1,[]);
y1_shape=size(y1);y1=reshape(y1,1,[]);
y2_shape=size(y2);y2=reshape(y2,1,[]);
y3_shape=size(y3);y3=reshape(y3,1,[]);
d_shape=size(d);d=reshape(d,1,[]);
w1_shape=size(w1);w1=reshape(w1,1,[]);
w2_shape=size(w2);w2=reshape(w2,1,[]);
w3_shape=size(w3);w3=reshape(w3,1,[]);
if isempty(k1p1), k1p1=0; end;
if isempty(k2p1), k2p1=0; end;
if isempty(k3p1), k3p1=0; end;
if isempty(k4p1), k4p1=0; end;
%
%***FIRST EXECUTABLE STATEMENT  TRI3
mm1 = fix(m - 1);
k1 = fix(k(1));
k2 = fix(k(2));
k3 = fix(k(3));
k4 = fix(k(4));
k1p1 = fix(k1 + 1);
k2p1 = fix(k2 + 1);
k3p1 = fix(k3 + 1);
k4p1 = fix(k4 + 1);
k2k3k4 = fix(k2 + k3 + k4);
if( k2k3k4~=0 )
l1 =fix(fix((k1+1)./(k2+1)));
l2 =fix(fix((k1+1)./(k3+1)));
l3 =fix(fix((k1+1)./(k4+1)));
lint1 = 1;
lint2 = 1;
lint3 = 1;
kint1 = fix(k1);
kint2 = fix(kint1 + k2);
kint3 = fix(kint2 + k3);
end;
for n = 1 : k1;
x = tcos(n);
if( k2k3k4~=0 )
if( n==l1 )
for i = 1 : m;
w1(i) = y1(i);
end; i = fix(m+1);
end;
if( n==l2 )
for i = 1 : m;
w2(i) = y2(i);
end; i = fix(m+1);
end;
if( n==l3 )
for i = 1 : m;
w3(i) = y3(i);
end; i = fix(m+1);
end;
end;
z = 1../(b(1)-x);
d(1) = c(1).*z;
y1(1) = y1(1).*z;
y2(1) = y2(1).*z;
y3(1) = y3(1).*z;
for i = 2 : m;
z = 1../(b(i)-x-a(i).*d(i-1));
d(i) = c(i).*z;
y1(i) =(y1(i)-a(i).*y1(i-1)).*z;
y2(i) =(y2(i)-a(i).*y2(i-1)).*z;
y3(i) =(y3(i)-a(i).*y3(i-1)).*z;
end; i = fix(m+1);
for ip = 1 : mm1;
i = fix(m - ip);
y1(i) = y1(i) - d(i).*y1(i+1);
y2(i) = y2(i) - d(i).*y2(i+1);
y3(i) = y3(i) - d(i).*y3(i+1);
end; ip = fix(mm1+1);
if( k2k3k4~=0 )
if( n==l1 )
i = fix(lint1 + kint1);
xx = x - tcos(i);
for i = 1 : m;
y1(i) = xx.*y1(i) + w1(i);
end; i = fix(m+1);
lint1 = fix(lint1 + 1);
l1 =fix(fix((lint1.*k1p1)./k2p1));
end;
if( n==l2 )
i = fix(lint2 + kint2);
xx = x - tcos(i);
for i = 1 : m;
y2(i) = xx.*y2(i) + w2(i);
end; i = fix(m+1);
lint2 = fix(lint2 + 1);
l2 =fix(fix((lint2.*k1p1)./k3p1));
end;
if( n==l3 )
i = fix(lint3 + kint3);
xx = x - tcos(i);
for i = 1 : m;
y3(i) = xx.*y3(i) + w3(i);
end; i = fix(m+1);
lint3 = fix(lint3 + 1);
l3 =fix(fix((lint3.*k1p1)./k4p1));
end;
end;
end; n = fix(k1+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;
tcos_shape=zeros(tcos_shape);tcos_shape(:)=tcos(1:numel(tcos_shape));tcos=tcos_shape;
y1_shape=zeros(y1_shape);y1_shape(:)=y1(1:numel(y1_shape));y1=y1_shape;
y2_shape=zeros(y2_shape);y2_shape(:)=y2(1:numel(y2_shape));y2=y2_shape;
y3_shape=zeros(y3_shape);y3_shape(:)=y3(1:numel(y3_shape));y3=y3_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
w1_shape=zeros(w1_shape);w1_shape(:)=w1(1:numel(w1_shape));w1=w1_shape;
w2_shape=zeros(w2_shape);w2_shape(:)=w2(1:numel(w2_shape));w2=w2_shape;
w3_shape=zeros(w3_shape);w3_shape(:)=w3(1:numel(w3_shape));w3=w3_shape;
end
%DECK TRIDIB