2D diffusion equation for a laser heat source

I am writing a code for laser source. However, when I am incorporating the source term, I am getting an error saying mismatch of array dimensions. I also want a plot between temperature and time. How to do that? Is there any issue with the CFL condition?

clc;
clear all;

% Parameters
rho=2600;
cp=473;

% domain and step
Lx=10;
Ly=10;

Nx=51;
Nt=3000;
Ny=51;
dx=Lx/(Nx-1);
dy=Ly/(Ny-1);

%CFL condition
c=1;
C=0.05;
dt=(C*dx/c);

%field variables
Tn=zeros(Ny,Nx);
x=linspace(0,Lx,Nx);
y=linspace(0,Ly,Ny);
[X,Y]=meshgrid(x,y);
q_dot=2.36*10^14*exp(-(800)*((x).^2+(y).^2));

K=ones(Ny,Nx);
K(1:Lx,1:Ly)=0.078;

%Initial condition
Tn(:,:)=0;
t=0;

%Loop
for n=1:Nt
    
    Tc=Tn;

    t=t+dt;
    for i=2:Nx-1
        for j=2:Ny-1
            Tn(j,i)=Tc(j,i)+ ...
            dt* (K(j,i)/rho/cp)*...
            ((Tc(j,i+1)+Tc(j+1,i)-4*Tc(j,i)+Tc(j-1,i)+Tc(j,i-1))/dx/dx);
        end
    end

    %Source term
    Sx=round(5.1*Nx/Lx);
    Sy=round(9.9*Ny/Ly);

    if (t<3)
        Tn(Sy,Sx)=Tn(Sy,Sx)+dt*q_dot/rho/cp;
    end

    %Boundary condition
    Tn(1,:)=0;
    Tn(end,:)=0;
    Tn(:,1)=0;
    Tn(:,end)=Tn(:,end-1);

    %Visualize
    mesh(X,Y,Tn);
    axis([0 Lx 0 Ly 0 50]);
    title(sprintf('Time=%f seconds',t))
    pause(0.1);
end

Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-51.

Answer 1

[Torsten]

qdot is 1x51, T(Sy,Sx) is a scalar. Thus

Tn(Sy,Sx)=Tn(Sy,Sx)+dt*q_dot/rho/cp;

tries to assign an array to a scalar which throws an error.

And why do you mix lengths and number of elements and switch between x and y in the definition of K ?

K=ones(Ny,Nx);
K(1:Lx,1:Ly)=0.078;

[AD]

Ik why I am getting the error? But how to fix it?

[Torsten]

In what region of your domain do you want to set the source term ? And how should the source term look like in this subdomain ?

"qdot" in your code is a vector of size 1x51, and you are lucky that Nx = Ny, otherwise MATLAB would have thrown an error when you compute "qdot" as

q_dot=2.36*10^14*exp(-(800)*((x).^2+(y).^2));

It is not obvious how "qdot" of size 1x51 is related to the domain which has a 51x51 grid.

Maybe you want

q_dot=2.36*10^14*exp(-(800)*((X).^2+(Y).^2));

I don't know.