newton raphson for nonlinear heat equation scheme

I am trying to apply newton raphson method to linearize a scheme arise from non linear heat equation of the form

and the scheme is :

where are non linear term that was splitted to convex-concave function.

I have written this code but I couldn't get what is my mistake and how to make it work.

I appreciate any one can help me with that

close all
clear all
clc 
%%
%PR parameters
R=8.314;
T=352.49;
% Critical constants for normal butane (c_6)
Tc = 507.44; % Critical temperature in Kelvin
Pc = 30.31; % Critical pressure in bar
w = 0.299;  % Acentric factor
% MW=0.05812;
m=0.37464+1.54226*w-0.26992*w^2;
a=0.45724*R^2*Tc^2/Pc*(1+m*(1-sqrt(T/Tc)))^2;
b=0.07780*R*Tc/Pc;
%Influenceparameter
m1=10^(-16)/(1.2326+1.3757*w);
m2=10^(-16)/(0.9051+1.5410*w);
c=a*b^(2/3)*(m1*(1-T/Tc)+m2);
%% Parameters
L = 0.000000008;        % Length of the domain
T = 0.1;      % Total time
Nx = 5;      % Number of spatial grid points
Nt = 10;    % Number of time steps
% c = 1;        % Thermal diffusivity constant
dx = L / (Nx - 1);  % Spatial step size
dt = T / Nt;  % Time step size

% Initial and boundary conditions
u0 = @(x) 7000-(7000/L)*x;  % Initial condition
uL = 7000;               % Dirichlet boundary condition at x=0
uR = 0;               % Dirichlet boundary condition at x=L

% Nonlinear source term and its convex/concave splitting
f = @(u) R.*T*(log(u./(1-u.*b))+u.*b./(1-u.*b))-u.*a./(1+2.*u.*b-u.^2*b.^2)-a./(2.*sqrt(2).*b).*log((1+(1-sqrt(2)).*u.*b)./(1+(1+sqrt(2)).*u.*b));        % Example: cubic nonlinearity
fc = @(u) R.*T.*(log(u./(1-u.*b))+u.*b./(1-u.*b));    % Convex part
fd = @(u) f(u) - fc(u);  % Concave part
fplot(f)

% Derivative of the convex part for Newton's method
dfc = @(u) R.*T*(1/(u.*(b.*u-1)).^2);   % Derivative of the convex part

% Exact solution function (if available)
% For demonstration, assuming an exact solution function
% exact_solution = @(x, t) 1./(1+exp(x-t)); 

% Spatial grid
x = linspace(0, L, Nx);
% u=zeros(Nt,Nx);
% Initial condition
u = u0(x);
figure(1)
plot(x,u0(x))
grid on
title('Initial condition')
xlabel('x')
ylabel('u0')

%%
% Laplacian matrix (central difference)
e = ones(Nx,1);
A = spdiags([e -2*e e], -1:1, Nx, Nx) / dx^2;

% Apply boundary conditions to A
A(1,:) = 7000; A(1,1) = 7000;   % Fix the boundary at x=0
A(end,:) = 0; A(end,end) = 0; % Fix the boundary at x=L 0;
%%
% Time integration loop
for n = 1:Nt
    % Apply boundary conditions at each time step
    u(1) = 0;   % Boundary at x=0
    u(Nx) = 0;  % Boundary at x=L

    % Initial guess for Newton's method (use the previous time step)
    U = u;
    
    % Newton iteration loop
    for k = 1:100
        % Calculate the residual F(U^k)
        F_Uk = U - dt * c * A * U - dt * fc(U) - (u + dt * fd(u));
        
        % Calculate the Jacobian J(U^k)
        J_Uk = speye(Nx) - dt * c * A - dt * spdiags(dfc(U), 0, Nx, Nx);
        
        % Solve for Newton update dU
        dU = J_Uk \ (-F_Uk);
        
        % Update U with the Newton step
        U = U + dU;
        
        % Check for convergence
        if norm(dU, inf) < 1e-6
            break;
        end
    end
    
    % Update solution for the next time step
    u = U;
    
    % Optionally, plot the solution at selected times
    if mod(n, Nt/10) == 0  % Plot every 10% of the total time steps
        figure;
        plot(x, u, '-o', 'DisplayName', 'Numerical Solution');
        xlabel('x');
        ylabel('u');
        title(sprintf('Solution at time t = %.2f', n * dt));
        grid on;
    end
end

Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.

Answer 1

[Torsten]

Why don't you use "pdepe" ? Ten lines of code maximum, and you are done.

To get rid of the error in the above code, replace

% Spatial grid
x = linspace(0, L, Nx);

by

% Spatial grid
x = linspace(0, L, Nx).';

But other problems follow.

Why do you set A as below ? Shouldn't this be u or U ?

A(1,:) = 7000; A(1,1) = 7000;   % Fix the boundary at x=0
A(end,:) = 0; A(end,end) = 0; % Fix the boundary at x=L 0;

Are you sure about u(1) = 0 ? I thought you want to keep it at 7000.

    % Apply boundary conditions at each time step
    u(1) = 0;   % Boundary at x=0
    u(Nx) = 0;  % Boundary at x=L

[Maryam Albuhayr]

yes you are right I want u(1)=7000