I am using spectral method to solve an ordinary differential equation. The following R code and the MATLAB code are completely the same. But the R code gets explode (the u value it solved would become very large), while MATLAB could stably solve the u. I do suspect this is because the computation stability of R is weaker than MATLAB. I would appreciate it if someone could analyze the situation in more details.

# R (The result would explode)
library(pracma)

f<- function(x){
  return(exp(x)*(x^2+ 4*x+1)+ exp(5*x)*(x^2-1)^5)
}
Cheb<- function(N){ #D_cheb of size N+1, corresponding to N sub-intervals
D= matrix(0, nrow= N+1, ncol=N+1)
x= cos((1:N)*pi/N)
for (i in (1:(N-1))){
  D[i+1,i+1]= -x[i]/(2*(1- x[i]^2))
}  
for(i in 1:(N-1)){
  for(j in 1:(N-1)){
    if(i !=j){
      D[i+1,j+1]= (-1)^(i+j)/(x[i]- x[j])
    }
  }
}
for (j in 2:N){
  D[1,j]= 2*(-1)^(j-1)/(1- x[j-1])
  D[N+1,j]= -2*(-1)^(j+N-1)/(1+ x[j-1])
  D[j,1]= -0.5*(-1)^(j-1)/(1- x[j-1])
  D[j,N+1]= 0.5*(-1)^(N+j-1)/(1+ x[j-1])
}
  D[1,1]= (2*N*N+1)/6; D[1,N+1]= 0.5*(-1)^N
  D[N+1,1]= -0.5*(-1)^N; D[N+1,N+1]= -(2*N*N+1)/6
  

return(D)

}

N= 64
D= Cheb(N)
A= D%*%D; A= A[2:N, 2:N]
x= cos((1:(N-1))*pi/N)
u= (x+1)^2*(x-1)
ff= vector(mode = "numeric", length = N-1)
for(i in 1:(N-1)){
  ff[i]= f(x[i])
}


while(max(abs(A%*%u+ u^5- ff))> 0.005){
  deltau= -inv(A+ diag(5*(u^4)))%*%(A%*%u+ u^5- ff)
  u= u+deltau
  print(max(abs(A%*%u+ u^5- ff)))
  print(u)
}

% MATLAB

f = @(x) exp(x).*(x.^2 + 4.*x + 1) + exp(5.*x).*(x.^2 - 1).^5;

function D = Cheb(N)
    % D_cheb of size N+1, corresponding to N sub-intervals
    D = zeros(N+1, N+1);
    x = cos((1:N)*pi/N);

    for i = 1:(N-1)
        D(i+1, i+1) = -x(i) / (2*(1 - x(i)^2));
    end

    for i = 1:(N-1)
        for j = 1:(N-1)
            if i ~= j
                D(i+1, j+1) = (-1)^(i+j) / (x(i) - x(j));
            end
        end
    end

    for j = 2:N
        D(1,   j)   =  2*(-1)^(j-1)      / (1 - x(j-1));
        D(N+1, j)   = -2*(-1)^(j+N-1)    / (1 + x(j-1));
        D(j,   1)   = -0.5*(-1)^(j-1)    / (1 - x(j-1));
        D(j,   N+1) =  0.5*(-1)^(N+j-1)  / (1 + x(j-1));
    end

    D(1,   1)   =  (2*N*N + 1)/6;
    D(1,   N+1) =  0.5*(-1)^N;
    D(N+1, 1)   = -0.5*(-1)^N;
    D(N+1, N+1) = -(2*N*N + 1)/6;
end

N  = 64;
D  = Cheb(N);
A  = D*D;
A  = A(2:N, 2:N);                  

x  = cos((1:(N-1))*pi/N).';        
u  = (x + 1).^2 .* (x - 1);        
ff = zeros(N-1, 1);

for i = 1:(N-1)
    ff(i) = f(x(i));
end

while max(abs(A*u + u.^5 - ff)) > 0.005
    deltau = -inv(A + diag(5*(u.^4))) * (A*u + u.^5 - ff);
    u= u + deltau;
end