%
% Demo for stablity 
%
clear all; close all;

N = 12;
%
% 解法一：向前递推，误差快速累积
S_f(1) = 0.182;  % S0
for n = 2 : N
    S_f(n) = 1/(n-1) - 5*S_f(n-1);
    S_f(n) = vpa(S_f(n),3);
end
disp('解法一：')
disp(S_f)

% 解法二：向后递推，误差得到很好的控制
S_b(N) = (1/(6*N) + 1/(5*N))/2;
S_b(N) = vpa(S_b(N),3);
for n = N-1 : -1 : 1
    S_b(n) = (1/n - S_b(n+1))/5;
    S_b(n) = vpa(S_b(n),3);
end
disp('解法二：')
disp(S_b)

% 精确值, 符号计算
syms x
for n = 1 : N
    St(n) = int(x^(n-1)/(x+5),0,1);
    St(n) = vpa(St(n),5);
end
disp('精确值：')
St = double(St);
disp(St)

%%% plot the figure 
plot(1:N,S_f,'vr', 1:N,S_b,'.b', 1:N,St,'ok');
legend('Foreward', 'Backward', 'Exact Sol');
