nlide_2

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

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

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

% Points for ploting.
points = linspace(x, 200); %#ok<NASGU>

% Looping for iterations.
for k = 1:5

% Derivate y1 and y2 using N matrix.
y1p = N*y1; y2p = N*y2;

% Compute right hand sides.
rhs1 = interporth(x, '1') + 0.5*chebypolyprod(y2p, y2p, y.n)';
rhs2 = interporth(x, '2*x');
rhs3 = -2*chebypolyprod(y1, y1p, y.n)';
rhs4 = -2*chebypolyprod(y2, y2p, y.n)';
rhs5 = -chebypolyprod(y1, y2p, y.n)'-chebypolyprod(y1p, y2, y.n)';

% Create the blocks for the T matrix. ---------------------------------------------------------
T11 = diff(y);
T12 = orthoval(x, M, 'coef', y2p)*diff(y) - volt(y, 'x-t');
T13 = zeros(y.n);
T14 = zeros(y.n);
T15 = volt(y);

T21 = volt(y, 'x-t');
T22 = diff(y);
T23 = volt(y);
T24 = volt(y);
T25 = zeros(y.n);

T31 = -2*(orthoval(x, M, 'coef', y1p)*y + orthoval(x, M, 'coef', y1)*diff(y));
T32 = zeros(y.n);
T33 = diff(y);
T34 = zeros(y.n);
T35 = zeros(y.n);

T41 = zeros(y.n);
T42 = -2*(orthoval(x, M, 'coef', y2p)*y + orthoval(x, M, 'coef', y2)*diff(y));
T43 = zeros(y.n);
T44 = diff(y);
T45 = zeros(y.n);

T51 = -orthoval(x, M, 'coef', y2p)*y - orthoval(x, M, 'coef', y2)*diff(y);
T52 = -orthoval(x, M, 'coef', y1p)*y - orthoval(x, M, 'coef', y1)*diff(y);
T53 = zeros(y.n);
T54 = zeros(y.n);
T55 = diff(y);
% ---------------------------------------------------------------------------------------------

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

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

% Truncating the blocks.
T(6: y.n+5-1, 1:y.n) = T11.mat(1:end-1, :);
T(6: y.n+5-1, y.n+1:2*y.n) = T12.mat(1:end-1, :);
T(6: y.n+5-1, 2*y.n+1:3*y.n) = T13(1:end-1, :);
T(6: y.n+5-1, 3*y.n+1:4*y.n) = T14(1:end-1, :);
T(6: y.n+5-1, 4*y.n+1:5*y.n) = -T15.mat(1:end-1, :);

T(y.n+5: 2*y.n+5-2, 1:y.n) = -T21.mat(1:end-1, :);
T(y.n+5: 2*y.n+5-2, y.n+1:2*y.n) = T22.mat(1:end-1, :);
T(y.n+5: 2*y.n+5-2, 2*y.n+1:3*y.n) = -T23.mat(1:end-1, :);
T(y.n+5: 2*y.n+5-2, 3*y.n+1:4*y.n) = T24.mat(1:end-1, :);
T(y.n+5: 2*y.n+5-2, 4*y.n+1:5*y.n) = T25(1:end-1, :);

T(2*y.n+5-1: 3*y.n+5-3, 1:y.n) = T31.mat(1:end-1, :);
T(2*y.n+5-1: 3*y.n+5-3, y.n+1:2*y.n) = T32(1:end-1, :);
T(2*y.n+5-1: 3*y.n+5-3, 2*y.n+1:3*y.n) = T33.mat(1:end-1, :);
T(2*y.n+5-1: 3*y.n+5-3, 3*y.n+1:4*y.n) = T34(1:end-1, :);
T(2*y.n+5-1: 3*y.n+5-3, 4*y.n+1:5*y.n) = T35(1:end-1, :);

T(3*y.n+5-2: 4*y.n+5-4, 1:y.n) = T41(1:end-1, :);
T(3*y.n+5-2: 4*y.n+5-4, y.n+1:2*y.n) = T42.mat(1:end-1, :);
T(3*y.n+5-2: 4*y.n+5-4, 2*y.n+1:3*y.n) = T43(1:end-1, :);
T(3*y.n+5-2: 4*y.n+5-4, 3*y.n+1:4*y.n) = T44.mat(1:end-1, :);
T(3*y.n+5-2: 4*y.n+5-4, 4*y.n+1:5*y.n) = T45(1:end-1, :);

T(4*y.n+5-3: 5*y.n, 1:y.n) = T51.mat(1:end-1, :);
T(4*y.n+5-3: 5*y.n, y.n+1:2*y.n) = T52.mat(1:end-1, :);
T(4*y.n+5-3: 5*y.n, 2*y.n+1:3*y.n) = T53(1:end-1, :);
T(4*y.n+5-3: 5*y.n, 3*y.n+1:4*y.n) = T54(1:end-1, :);
T(4*y.n+5-3: 5*y.n, 4*y.n+1:5*y.n) = T55.mat(1:end-1, :);

b(0*y.n+5+1: 1*y.n+4) = rhs1(1:end-1);
b(1*y.n+5-0: 2*y.n+3) = rhs2(1:end-1);
b(2*y.n+5-1: 3*y.n+2) = rhs3(1:end-1);
b(3*y.n+5-2: 4*y.n+1) = rhs4(1:end-1);
b(4*y.n+5-3: 5*y.n) = rhs5(1:end-1);

% Solve the linear system.
a = T\b;  a = reshape(a, y.n, 5);
y1 = a(:, 1);  y2 = a(:, 2);  y3 = a(:, 3);  y4 = a(:, 4);  y5 = a(:, 5);

%figure(1)
%subplot(1, 2, 1)
%semilogy(points, abs(orthoval(x, points, 'coef', y1) - (sinh(points))));
%hold on
%subplot(1, 2, 2)
%semilogy(points, abs(orthoval(x, points, 'coef', y2) - (cosh(points))));
%hold on

end

figure(2)
 subplot(1, 2, 1)
 xlabel('$x$', 'Interpreter', 'latex')
 ylabel('$|\sinh(x)-y_1^{(k)}(x)|$', 'Interpreter', 'latex')
 text(0.5, 1e-10, 'k = 1', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')
 text(0.5, 1e-10, 'k = 2', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')
 text(0.5, 1e-10, 'k = 3', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')
 text(0.5, 1e-10, 'k = 4', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')
 text(0.5, 1e-10, 'k>4', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')

 subplot(1, 2, 2)
 xlabel('$x$', 'Interpreter', 'latex')
 ylabel('$|\cosh(x)-y_2^{(k)}(x)|$', 'Interpreter', 'latex')
 text(0.5, 1e-10, 'k = 1', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')
 text(0.5, 1e-10, 'k = 2', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')
 text(0.5, 1e-10, 'k = 3', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')
 text(0.5, 1e-10, 'k>3', 'Backgroundcolor', [1 1 1], 'Edgecolor', 'k')

data = [0 0 0
0.1 1.39e-17 0
0.2 5.27e-16 0
0.3 4.45e-14 1.11e-15
0.4 1.05e-12 3.51e-14
0.5 1.23e-11 5.10e-13
0.6 9.11e-11 4.55e-12
0.7 4.97e-10 2.90e-11
0.8 2.16e-9  1.44e-10
0.9 7.90e-9  5.92e-10
1   2.52e-8  2.10e-9];

points = 0:0.01:1;
data4 = abs(orthovalv(y1, points, x.domain, x.basis) - (sinh(points)));
data5 = abs(orthovalv(y2, points, x.domain, x.basis) - (cosh(points)));

figure(2)
subplot(1, 2, 1)
semilogy(points, abs(orthovalv(y1, points, x.domain, x.basis) - (sinh(points))), ...
         'LineWidth', 1.6); hold on
plot(data(:, 1), data(:, 2), 'LineWidth', 1.6), hold on
xlabel('$x$', 'Interpreter', 'latex', 'FontSize', 16)
ylabel(['$|\sinh(x)-y_1^{(',num2str(k),')}(x)|$'], 'Interpreter', 'latex', 'FontSize', 16)
legend('Tau Toolbox', 'AbTa09')
subplot(1, 2, 2)
semilogy(points, abs(orthovalv(y2, points, x.domain, x.basis) - (cosh(points))), ...
         'LineWidth', 1.6); hold on
plot(data(:, 1), data(:, 3), 'LineWidth', 1.6), hold on
xlabel('$x$', 'Interpreter', 'latex', 'FontSize', 16)
ylabel(['$|\cosh(x)-y_2^{(',num2str(k),')}(x)|$'], 'Interpreter', 'latex', 'FontSize', 16)
legend('Tau Toolbox', 'AbTa09')