crank nicolson scheme for finite difference method in matlab code

Dear sir/madam

I am trying to solve the finite difference methof for crank nicolson scheme to 2d heat equation. please let me know if you have any MATLAB CODE for this

Boundary codition are

If you can kindly send me the matlab code, it will be very useful for my research work . thank you very much.

% Equation of energy

clear; clc; clear all;

rho=1;

cp=1;

Lx=8;

Ly=1;

Nx=26;

%Nt=5;

Ny=26;

dx=0.05;%%Lx/(Nx-1);

dy=0.05;%%%Ly/(Ny-1);

c=1;

C=0.05;

dt=0.01;%%%C*dx/c;

Tk=zeros(Nx,Ny);

Uk=zeros(Nx,Ny);

Ck=zeros(Nx,Ny);

Vk=zeros(Nx,Ny);

x=linspace(0,Lx,Nx);

y=linspace(0,Ly,Ny);

[x,y]=meshgrid(x,y);

M=0.5;

N=0.4;

Pr=0.71;

eps=1.0;

del=1;

Uk(:,:)=0;

Tk(:,:)=20;

Ck(:,:)=0;

Vk(:,:)=0;

t=0;

tol = 1e-6;

error = 1;

k = 0;

while error > tol

k = k+1;

Tkold = Tk;

for i = 2:Nx-1;

for j = 2:Ny-1;

%Equation of Energy

Tk(i,j)= (Tk(i,j-1)*[-Vk(i,j)*dt/4*dy-1/Pr*(3*N+4/3*N)*dt/2*dy*dy])...

-Tk(i,j)*[1+Uk(i,j)*dt/2*dx+1/Pr*((3*N+4)/3*N)*dt/dy*dy-del*(dt/2)]...

-Tk(i,j+1)*[Vk(i,j)*dt/4*dy-1/Pr*(3*N+4/3*N)*dt/2*dy*dy]...

-Uk(i,j)*[(Tk(i-1,j)+Tk(i-1,j)-Tk(i,j)*dt)/2*dx]...

+Vk(i,j)*[(Tk(i,j-1)+Tk(i,j+1)*dt)/4*dy]...

+(1/Pr)*((3*N+4)/3*N)*[(Tk(i,j-1)-2*Tk(i,j)+Tk(i,j+1))*dt/2*dy*dy]...

+del*Tk(i,j)*(dt/2)-eps*[(Uk(i,j+1)-Uk(i,j))/dy].^2;

end

end

error = min(min(abs(Tkold-Tk)));

end

datavalues=[x' y' Uk' Vk' Ck' Tk'];

figure(1)

subplot(3,1,1), contour(x,y,Tk), colormap

title('Temperature (transient state)'), xlabel('x'), ylabel('y'),colorbar

subplot(3,1,2), pcolor(x,y,Tk), shading interp,

title('Temperature (transient state)'), xlabel('x'), ylabel('y'),colorbar

subplot(3,1,3)

surf(Tk')

xlabel('x')

ylabel('y')

zlabel('U')

colorbar

figure(2)

plot(y,Tk)

figure(3)

subplot(3,1,1), contour(x,y,Tk), colormap

title('Temperature (steady state)'), xlabel('x'), ylabel('y'),colorbar

subplot(3,1,2), pcolor(x,y,Tk), shading interp,

title('Temperature (steady state)'), xlabel('x'), ylabel('y'),colorbar

subplot(3,1,3)

surf(Vk')

xlabel('x')

ylabel('y')

zlabel('V')

colorbar

kindly correct the matlab code. thank you,