Implicit method in slab

Hi every one, This is my first code I have ever written with some help. I just would like to make sure is good and right before I submitted to my instructor.

The assignment is to solve a heat diffusion in a slab with two neumann boundaries (convection) and one initial (30 C).

I have attached a file for the main equations.

Thanks in advance

% [u,err,x,t] = heat2(t_0,t_f,M,N)

% function [u,err,x,t] = heat2(t_0,t_f,M,N)

% define the mesh in space

N=101; %Number of space intervals M=2000; %number of time steps l=1; %thickness of wall in meter;; t_0=0; %initial time t_f=10; %final time (sec) dx = l/(N-1); %size of the mesh x = 0:dx:N-1; x = x'; %inverse of mesh

bi=0.77; %coefficient of terms bo=0.88;

Te=34; %exterior wall temperature k=0.7; %thermal conductivity of bricks roh=1600; %density of bricks (kg/m^3) Cp=840; %Cp of bricks (J/kg-K) T_in=22; %Inner wall temperature alfa=Cp*roh/k;

% define the mesh in time dt = (t_f-t_0)/M; t = t_0:dt:t_f;

% define the ratio r r=alfa*dt/dx^2; %Fourier number

Beta=[2*r*bo*Te; zeros(N-2,1); 2*r*bi*T_in];

%initial condition for i=1:N %u(i,1) = cos(x(i)); u(i,1) = 15; end err(:,1) = u(:,1) - exp(t_0-t(1))*cos(x);

A = zeros(N,N); A(1,1) = 1+2*r+2*r*bo; A(1,2) = -2*r; % A(1,3:N) = 0; for i=2:N-1 A(i,i-1) = -r; A(i,i) = 1+2*r; A(i,i+1) = -r; end %A(N,N+1) = 0; A(N,N-1) = -2*r; A(N,N) = 1+2*r+2*r*bi;

[L,UU] = LU(A);

%[y] = down_solve(L,u_old); %[u_new] = up_solve(UU,y);

for j=1:M %time step [y] = down_solve(L,u(:,j)+Beta); u(:,j+1) = up_solve(UU,y); err(:,j+1) = u(:,j+1) - exp(t_0-t(j+1))*cos(x);

the sub_files:

1)

% [L,U] = LU(A) LU

function [L,U] = LU(A);

[m,n] = size(A);

for i=1:n L(i,i) = 1; end

for k=1:n-1 for i=k+1:n mult = A(i,k)/A(k,k); for j=k+1:n A(i,j) = A(i,j) - A(k,j)*mult; end L(i,k) = mult; end end

for i=1:n for j=i:n U(i,j) = A(i,j); end end

2) % [x] = up_solve(U,b) takes an upper triangular matrix U and vector b % and solves Ux=b

function [x] = up_solve(U,b)

[m,n] = size(U);

x(n) = b(n)/U(n,n);

for k=1:n-1 x(n-k) = b(n-k); for j=n-k+1:n x(n-k) = x(n-k) - U(n-k,j)*x(j); end x(n-k) = x(n-k)/U(n-k,n-k); end x = x';

3)

% [x] = down_solve(L,b) takes a lower triangular matrix L and vector b % and solves Lx=b

function [x] = down_solve(L,b)

[m,n] = size(L);

x(1) = b(1)/L(1,1);

for k=2:n x(k) = b(k); for j=1:k-1 x(k) = x(k) - L(k,j)*x(j); end x(k) = x(k)/L(k,k); end x = x';