syms x k;
N=25; %order of the polynomial, N>1
%1st interval [1,e], center of U at 1
U=zeros(1,N+1); %inicialization of U
U(1)=exp(1); %initial value
F=@(k)((-1).^(k+1))./k; %DT of ln at 1
FF=[0 F(1:N-1)]; %values of F
E=@(k) 1./factorial(k); %DT of exp at 0
EE=[1 E(1:N-1)]; %values of E
eqn1='(k+1)*x==U(k+1)-kd(k-1)-kd(k)+H'; % recurrent equation, x at the place of U(k+2)
B=zeros(N,N); %inicialization of a matrix with Bell's coefficients 
B(1,1)=1; %first value of B 
k=0; %first step (no need to use B)
H=1; %H[1]=E[0]
U(k+2)=solve(eval(eqn1),x); %U(2)
for k=1:N-1 %calculation of values U(3) - U(N+1)
  for l=2:k+1 % fill (k+1)th row of Bell's matrix in columns 2 to diagonal
  for i= 1:k+2-l
  B(k+1,l)=B(k+1,l)+((i*(l-1))/(k))*FF(i+1)*B(k+1-i,l-1);
  end
  end
  H=dot(EE,B(k+1,:)); %H(k) scalar product of values of F and (k+1)th row of matrix B
    U(k+2)=solve(eval(eqn1),x); %get next member of vector U
end
sol1=dot(U,(x-1).^(0:N)); %solution in the form of polynomial in x