set up coefficient matrix and rhs vector for a 2d unsteady heat conduction equation. use implicit euler in time and central differencing in space. then use jacobi solver in matlab

main.m
% set up parameters
L = 1; % length of domain
H = 1; % height of domain
T = 1; % total time
dx = 0.05; % grid spacing in x direction
dy = dx; % grid spacing in y direction
dt = 0.0001; % time step
alpha = 0.1; % thermal diffusivity

% set up grid
x = 0:dx:L;
y = 0:dy:H;
[X,Y] = meshgrid(x,y);

% set initial condition
u0 = zeros(length(x),length(y));
u0(1,:) = 1; % boundary condition

% set up coefficient matrix
nx = length(x);
ny = length(y);
A = zeros(nx*ny);
for i = 2:nx-1
    for j = 2:ny-1
        n = j + (i-1)*ny;
        A(n,n) = 1 + 4*alpha*dt/dx^2 + 4*alpha*dt/dy^2;
        A(n,n-1) = -alpha*dt/dx^2;
        A(n,n+1) = -alpha*dt/dx^2;
        A(n,n-ny) = -alpha*dt/dy^2;
        A(n,n+ny) = -alpha*dt/dy^2;
    end
end

% set up right-hand side vector
F = zeros(nx*ny,1);

% solve matrix equation using Jacobi method
u = u0(:);
for t = dt:dt:T
    u_star = u;
    for k = 1:100 % maximum of 100 iterations
        for i = 2:nx-1
            for j = 2:ny-1
                n = j + (i-1)*ny;
                F(n) = u_star(n) + alpha*dt/dx^2*(u_star(n-1) - 2*u_star(n) + u_star(n+1)) ...
                    + alpha*dt/dy^2*(u_star(n-ny) - 2*u_star(n) + u_star(n+ny));
            end
        end
        u = A\F;
        if norm(u - u_star) < 1e-6 % convergence criterion
            break;
        end
    end
end

% reshape solution
u = reshape(u,nx,ny);

% plot solution
surf(X,Y,u)
xlabel('x')
ylabel('y')
zlabel('temperature')
1446 chars
64 lines