% script recur
% to find minimal solution of
% y(n+1) + a(n)*y(n) + b(n)*y(n-1) = ( f(n) ); 
% set y(0) = 1; y(N+K) = 0; K buffer where coefficients are unconverged
% or
%A=/ a(1)   1     0      0    .... \  d= / -b(1) \
%  | b(2)  a(2)   1      0    .... |     |   0   |
%  |  0    b(3)  a(3)    1    .....| ,   |   .   |
%  |  0     0    b(4)   a(4)  .....|     |   .   |
%  | .....                  .....  |     |   .   |
%  \ .....              b(N)   a(N)/     \   0   /
%  solve A*y = d
format long e;
x =1; N = 20; k = 5;
e = -2/x;
for i = 1:N+k
    a(i) = e*i;
    b(i) = 1;
    c(i) = 1;
end
% solution as tridiagonal system
A = diag(a)+diag(b(2:N+k),-1)+diag(c(1:N+k-1),1);
d = zeros(N+k,1);
d(1,1) = -b(1)*.7651976866;
y = A\d;
y = [.7651976866;y];
% solution with Miller's algorithm
r = zeros(N+k,1);
r(N+k+1,1) = 0;
for i = N+k-1:-1:0
    r(i+1,1) = -b(i+1)/(a(i+1)+r(i+2,1));
end
y1(1,1) = .7651976866;
prod = r(1,1);
for i = 2:N+k+1
    y1(i,1) = prod*y1(1,1);
    prod = prod*r(i,1);
end
y
y1