% Create tau objects.
[x, y] = tau('ChebyshevT', [0 1], 129);

% Vector with initial approximations (satisfing the conditions).
y1 = zeros(y.n, 1); y1(1) = 1;
y2 = y1;

% Matrix N and M in the Chebyshev basis.
N = matrixN(y);  M = matrixM(x);

% Compute right hand side 1.
rhs1 = interporth(x, '0.5*(exp(-2)-1)');

% Compute blocks independent of the iteration.
T11 = diff(y)+y;  T12 = fred(y);  T22 = diff(y);

% Memory allocation for T matrix ans b vector.
T = zeros(2*y.n);  b = zeros(2*y.n, 1);

% Write the conditions.
T(1, 1:y.n)       = orthoval(x, 0, 'difforder', 0); b(1) = 1; % y1(0) = 1
T(2, y.n+1:2*y.n) = orthoval(x, 0, 'difforder', 0); b(2) = 1; % y2(0) = 1

% Blocks truncation.
T(3: y.n+2-1, 1:y.n) = T11.mat(1:end-1, :);
T(3: y.n+2-1, y.n+1:2*y.n) = -T12.mat(1:end-1, :);
T(y.n+2: 2*y.n+2-2, y.n+1:2*y.n) = T22.mat(1:end-1, :);

% Allocate the right hand side 1.
b(3: y.n+2-1) = rhs1(1:end-1);

% Points for ploting.
points = linspace(x, 200);

% Loop for iterations.
for k = 1:6

% Derivate y1 using N matrix.
y1p = N*y1;

% Compute right hand side 2 and allocate space.
rhs2 = -2*chebypolyprod(y1, y1p, y.n)';
b(y.n+2: 2*y.n+2-2) = rhs2(1:end-1);

% Block truncation depending on iteration and allocate space.
T21 = -2*(orthoval(x, M, 'coef', y1p)*y + orthoval(x, M, 'coef', y1)*diff(y));
T(y.n+2: 2*y.n+2-2, 1:y.n) = T21.mat(1:end-1, :);

% Solve the linear system.
a = T\b;  a = reshape(a, y.n, 2); y1 = a(:, 1);  y2 = a(:, 2);

semilogy(points, abs(orthoval(x, points, 'coef', y1) - (exp(-points))), ...
                                                    'LineWidth', 1.6, 'color', colorR2B(6, k));
hold on
end

xlabel('$x$', 'Interpreter', 'latex')
ylabel('$|\exp(-x)-y^{(k)}(x)|$', 'Interpreter', 'latex')

Published with MATLAB® R2016b