%% Using Function Handle Instead of Numeric Matrix
% Solve a linear system by providing bicgstabl with a function handle that computes A*x in place of the coefficient matrix A.
%
% One of the Wilkinson test matrices generated by gallery is a 21-by-21 tridiagonal matrix. Preview the matrix.

A = gallery('wilk',21)

%
% The Wilkinson matrix has a special structure, so you can represent the operation A*x with a function handle. When A multiplies a vector, most of the elements in the resulting vector are zeros. The nonzero elements in the result correspond with the nonzero tridiagonal elements of A. Moreover, only the main diagonal has nonzeros that are not equal to 1.
% The expression \textrm{Ax} becomes:
% \textrm{Ax}=\left\lbrack \begin{array}{cccccccccc}
% 10 & 1 & 0 & \cdots  &  & \cdots  &  & \cdots  & 0 & 0\\
% 1 & 9 & 1 & 0 &  &  &  &  &  & 0\\
% 0 & 1 & 8 & 1 & 0 &  &  &  &  & \vdots \\
% \vdots  & 0 & 1 & 7 & 1 & 0 &  &  &  & \\
%  &  & 0 & 1 & 6 & 1 & 0 &  &  & \vdots \\
% \vdots  &  &  & 0 & 1 & 5 & 1 & 0 &  & \\
%  &  &  &  & 0 & 1 & 4 & 1 & 0 & \vdots \\
% \vdots  &  &  &  &  & 0 & 1 & 3 & \ddots  & 0\\
% 0 &  &  &  &  &  & 0 & \ddots  & \ddots  & 1\\
% 0 & 0 & \cdots  &  & \cdots  &  & \cdots  & 0 & 1 & 10
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_1 \\
% x_2 \\
% x_3 \\
% x_4 \\
% x_5 \\
% \vdots \\
% \\
% \vdots \\
% \\
% x_{21}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% 10x_1 +x_2 \\
% x_1 +9x_2 +x_3 \\
% x_2 +8x_3 +x_4 \\
% \vdots \\
% x_{19} +9x_{20} +x_{21} \\
% x_{20} +10x_{21}
% \end{array}\right\rbrack.
% The resulting vector can be written as the sum of three vectors:
% \textrm{Ax}=\left\lbrack \begin{array}{c}
% 0+10x_1 +x_2 \\
% x_1 +9x_2 +x_3 \\
% x_2 +8x_3 +x_4 \\
% \vdots \\
% x_{19} +9x_{20} +x_{21} \\
% x_{20} +10x_{21} +0
% \end{array}\right\rbrack=\left\lbrack \begin{array}{c}
% 0\\
% x_1 \\
% \vdots \\
% x_{20}
% \end{array}\right\rbrack +\left\lbrack \begin{array}{c}
% 10x_1 \\
% 9x_2 \\
% \vdots \\
% 10x_{21}
% \end{array}\right\rbrack +\left\lbrack \begin{array}{c}
% x_2 \\
% \vdots \\
% x_{21} \\
% 0
% \end{array}\right\rbrack.
% In MATLAB®, write a function that creates these vectors and adds them together, thus giving the value of A*x:
% function y = afun(x)
% y = [0; x(1:20)] + ...
%     [(10:-1:0)'; (1:10)'].*x + ...
%     [x(2:21); 0];
% end
% (This function is saved as a local function at the end of the example.)
%
% Now, solve the linear system \textrm{Ax}=b by providing bicgstabl with the function handle that calculates A*x. Use a tolerance of 1e-12 and 50 iterations.

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = bicgstabl(@afun,b,tol,maxit)

%
% Check that afun(x1) produces a vector of ones.

afun(x1)

%
%% Local Functions

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

% Copyright 2019 The MathWorks, Inc.
%