There are some mistakes in your article.
The code is based on
% Set parameters
X0 = 0.0;
Y0 = 1.0;
Z0 = 0.0;
N = 120;
SIGMA = 10;
R = 28;
B = 8/3;
% Initialize coefficient arrays
a = zeros(N+1,1);
b = zeros(N+1,1);
c = zeros(N+1,1);
a(1) = X0;
b(1) = Y0;
c(1) = Z0;
% Compute Taylor coefficients
for k = 1:N
SB = 0.0;
SC = 0.0;
for i= 0:k-1
SB = SB + a(i+1)*c(k-i);
SC = SC + a(i+1)*b(k-i);
end
a(k+1) = (SIGMA*(b(k) - a(k)))/k;
b(k+1) = ((R*a(k) - b(k)) - SB)/k ;
c(k+1) = (-B*c(k) + SC)/k ;
end
% Evaluate Taylor series at T
T = 0:0.01:0.2;
X = zeros(size(T));
Y = zeros(size(T));
Z = zeros(size(T));
for i = 1:numel(T)
tau = T(i);
x = a(N+1);
y = b(N+1);
z = c(N+1);
for k = N:-1:1
x = x*tau + a(k);
y = y*tau + b(k);
z = z*tau + c(k);
end
X(i) = x;
Y(i) = y;
Z(i) = z;
end
figure(1)
plot(T,X,'Color','red')
hold on
figure(2)
plot(T,Y,'Color','red')
hold on
figure(3)
plot(T,Z,'Color','red')
hold on
% Compare with ODE solution
fun = @(t,x)[SIGMA*(x(2)-x(1));x(1)*(R-x(3))-x(2);x(1)*x(2)-B*x(3)];
x0 = [X0;Y0;Z0];
tspan = T;
options = odeset ("RelTol", 1e-8, "AbsTol", 1e-8);
[TT,XX] = ode15s(fun,tspan,x0,options);
figure(1)
plot(TT,XX(:,1),'Color','blue')
figure(2)
plot(TT,XX(:,2),'Color','blue')
figure(3)
plot(TT,XX(:,3),'Color','blue')
