function [x]= penta2(y,a,b,c)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Kurt Annen annen@web-reg.de
% Date: 15/05/2004
% Internet: www.web-reg.de
%
% Solves the problem Ax=b when A is pentadiagonal and strongly nonsingular.
% This is much faster than x=A\y for large matrices.
%
% Reference: Sp䴨, Helmuth "Numerik: Eine Einf? f?hematiker und Informatiker"
% S. 110 . Vieweg-Verlag Braunschweig/Wiesbaden (1994)
%
% a = main diagonal
% b = 2. diagonal
% c = 3. diagonal
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin ~= 4
error('penta(A,y) requires four arguments');
end
[n] = length(a);
[m] = length(y);
[o] = length(b);
[p] = length(c);
if m ~= n
error('Inner matrix dimensions must agree');
end
c(m)=0;
c(m-1)=0;
b(m)=0;
if m ~= length(b) & m ~= length(c)
error('a,b,c must have the same dimension');
end
% a optimized algorithm
h1=0;
h2=0;
h3=0;
h4=0;
h5=0;
hh1=0;
hh2=0;
hh3=0;
hh4=0;
hh5=0;
z=0;
hb=0;
hc=0;
for i=1:m
z=a(i)-h4*h1-hh5*hh2;
hb=b(i);
hh1=h1;
h1=(hb-h4*h2)/z;
b(i)=h1;
hc=c(i);
hh2=h2;
h2=hc/z;
c(i)=h2;
a(i)=(y(i)-hh3*hh5-h3*h4)/z;
hh3=h3;
h3=a(i);
h4=hb-h5*hh1;
hh5=h5;
h5=hc;
end
h2=0;
h1=a(m);
y(m)=h1;
for i=m:-1:1
y(i)=a(i)-b(i)*h1-c(i)*h2;
h2=h1;
h1=y(i);
end
x=y;
end
return