eigen value of stochastic matrix

I wan to find eigen values and eigen vector of a stochastic matrix using Newton Raphson method.

% Define the dimensions and variables
n = 3; % Size of vectors and matrices
m = 2; % Number of matrices
P = 4; % Number of eigenvalues and tensors
maxIterations = 100; % Maximum number of iterations
tolerance = 1e-6; % Tolerance for convergence

% Initialize variables
u = ones(n*P, 1); % Initial guess for u, each column represents u1, u2, ..., uP
lambda = ones(P, 1); % Initial guess for lambda, each element represents lambda1, lambda2, ..., lambdaP

% Define the matrices and tensors
c0 = eye(P); % Example: identity matrix
A0 = magic(n); % Example: identity matrix
c = rand(P, P, m); % Example: random P x P x m tensor
A = rand(n, n, m); % Example: random n x n x m tensor
e = rand(P, P, P); % Example: random P x P x P tensor
cA=0;
for i=1:m
    cA= cA+kron(c(:,:,i),A(:,:,i));
end

% Perform Newton-Raphson iteration
for iter = 1:maxIterations
    % Compute the left-hand side of the equation
    
    lhs = (kron(c0 , A0) + cA) * u;
    
    % Compute the right-hand side of the equation
    
    lame=0;
    for p = 1:P
        
        lame = lame+kron(lambda(p) * e(:, :, p),eye(n)) ;
    end
    rhs=lame*u
    % Compute the residual and check for convergence
    residual_1 = lhs - rhs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize the variables
I=eye(n,n);
cA1=zeros(P,1);


    for i=1:P
    cA1(i)= u'*kron(e(:,:,i),I)*u;
    end

   
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
residual_2=cA1;

% Compute the residual and check for convergence
    residual = [residual_1;residual_2];
    if norm(residual(:)) < tolerance
        disp('Converged!');
        break;
    end
   
   jacobian = jac_1(A,lambda , u, n,m,P);
    
    % Update u and lambda using the Newton-Raphson method
    delta = jacobian \ residual(:);
    u = u - delta(1:n*P);
    
    lambda = lambda - delta(n*P+1:end);
    
    % Display the current iteration and residual
    disp(['Iteration: ' num2str(iter) ', Residual: ' num2str(norm(residual(:)))]);
end

uf=reshape(u,n,[])
lambdaf=lambda

but the answer I am getting is

lambdaf =

1.0e+58 *

2.7709

-1.9052

0.3483

-0.3253.

It would be very helpful if someone help me improve this code

[Kautuk Raj]

Can you tell us what answer are you expecting?

[rakesh kumar]

the lambdaf and uf values should be small.

tlength=1;
nel=2;


kmn=1;
w=0.95;
sigma=0.6;

E0=2.0*10^11;
el=E0*ones(1,nel);
[alpha,kl]=klterm4_spanos3(w,sigma,kmn,nel);

[kk,mm]=stiffm_new(nel,el);
kk=kk(3:end,3:end);
[ka,ma]=stiffm(nel,1);

m=ma(3:end,3:end);
for t=1:kmn
for j=1:kl
[kk2(:,:,j)]=stiffm_new(nel,E0*alpha(j,:,t));
kk1(:,:,j)=kk2(3:end,3:end,j);
end

end
k00=reshape(kk1,2*nel,[]);
k=[kk k00];
k=reshape(k,2*nel,2*nel,[]);

M=kl;

p=1;

A0=k(:,:,1);
A=k(:,:,2:end);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Define the dimensions and variables
n = 2*nel; % Size of vectors and matrices
%M = 2; % Number of matrices
%p = 2; % order of polynomials
P=factorial(M+p)/(factorial(M)*factorial(p));
m1 = 1; % m1th eigen value
maxIterations = 1; % Maximum number of iterations
tolerance = 1e-6; % Tolerance for convergence

% Initialize variables
u1 = ones(n*P, 1); % Initial guess for u, each column represents u1, u2, ..., uP
lambda1 = ones(P, 1); % Initial guess for lambda, each element represents lambda1, lambda2, ..., lambdaP

method=1;
c00= c_ijk(M,p,method);
e= e_ijk(M,p);
c_ijk_mat=my_cell2mat(c00,M,p);
e1=my_cell2mat_e(e,M,p);

% Define the matrices and tensors
%c0 = eye(P); % Example: identity matrix
c0 = c_ijk_mat(:,:,1);
%A0 = magic(n); % Example: identity matrix
%c = rand(P, P, M); % Example: random P x P x m tensor
c = c_ijk_mat(:,:,2:end);
%A = 10^-2*rand(n, n, M); % Example: random n x n x m tensor
% e = rand(P, P, P); % Example: random P x P x P tensor
e= e1;

[newLambda, newU] = replaceLambdaAndU(m1,P,A0);
lambda=newLambda;
u=reshape(newU,[],1);

cA=0;
for i=1:M
    cA= cA+kron(c(:,:,i),A(:,:,i));
end

% Perform Newton-Raphson iteration
for iter = 1:maxIterations
    % Compute the left-hand side of the equation
    
    lhs = (kron(c0 , A0) + cA) * u;
    
    % Compute the right-hand side of the equation
    
    lame=0;
    for p = 1:P
        
        lame = lame+kron(lambda(p) * e(:, :, p),eye(n)) ;
    end
    rhs=lame*u;
    % Compute the residual and check for convergence
    residual_1 = lhs - rhs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize the variables
I=eye(n,n);
cA1=zeros(P,1);


    for i=1:P
    cA1(i)= u'*kron(e(:,:,i),I)*u;
    end

   
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
residual_2=cA1;

% Compute the residual and check for convergence
    residual = [residual_1;residual_2];
    if norm(residual(:)) < tolerance
        disp('Converged!');
        break;
    end
   
   jacobian = jac_1(A,lambda , u, n,M,P);
    
    % Update u and lambda using the Newton-Raphson method
    delta = jacobian \ residual(:);
    u = u - delta(1:n*P);
    
    lambda = lambda - delta(n*P+1:end);
    
    % Display the current iteration and residual
    disp(['Iteration: ' num2str(iter) ', Residual: ' num2str(norm(residual(:)))]);
end

uf=reshape(u,n,[])
lambdaf=lambda

I want first elemnets of lambdaf vetor to be smallest value of eigen value of A0 and all other values of lambdaf to be smaller values than lambdaf ,of the order of 10^-2*lambdaf