Finish the code heat$2$Dimplicit.m $($shared in Files$),$ by programming the implicit finite difference approximation of the $2$D temperature equation.

$\%$ Solves the $2$D heat equation with an implicit finite difference scheme

clear

$\%$Physical parameters

L $=150$e$3$; $\%$ Width of lithosphere $[$m$]$

H $=100$e$3$; $\%$ Height of lithosphere $[$m$]$

Tbot $=1300$; $\%$ Temperature of bottom lithosphere $[$C$]$

Tsurf $=0$; $\%$ Temperature of country rock $[$C$]$

Tplume $=1500$; $\%$ Temperature of plume $[$C$]$

kappa $=1$e$-6$; $\%$ Thermal diffusivity of rock $[$m$2/$s$]$

Wplume $=25$e$3$; $\%$ Width of plume $[$m$]$

day $=3600*24$; $\%$ # seconds per day

year $=365\mathrm{.}25*$day; $\%$ # seconds per year

dt $=100$e$6*$year; $\%$ timestep

$\%$ Numerical parameters

nx $=51$; $\%$ # gridpoints in x$-$direction

nz $=51$; $\%$ # gridpoints in z$-$direction

nt $=100$; $\%$ Number of timesteps to compute

dx $=$ L$/($nx$-1)$; $\%$ Spacing of grid in x$-$direction

dz $=$ H$/($nz$-1)$; $\%$ Spacing of grid in z$-$direction

$[$x$2$d$,$z$2$d$]=$ meshgrid$(-$L$/2$:dx:L$/2,-$H:dz:$0)$; $\%$ create grid

$\%$ Setup initial linear temperature profile

T $=$ abs$($z$2$d$./$H$)*$Tbot;

$\%$ Imping plume beneath lithosphere

ind $=$ find$($abs$($x$2$d$(1,$:$))<=$ Wplume$/2)$;

T$(1,$ind$)=$ Tplume;

$\%$ Setup numbering

num $=1$;

for i$=1$:nz

for j$=1$:nx

Number$($i$,$j$)=$ num;

num $=$ num$+1$;

end

end

$\%$ Construct the A matrix

A $=$ sparse$($nx$*$nz$,$nx$*$nz$)$;

sx $=$ kappa$*$dt$/$dx$^2$;

sz $=$ kappa$*$dt$/$dz$^2$;

for i $=2$:nz$-1$

for j $=2$:nx$-1$

ii $=$ Number$($i$,$j$)$;

A$($ ii$,$ Number$($i$+1,$j $))=??$;

A$($ ii$,$ Number$($i $,$j$+1))=??$;

$??$

end

end

$\%$ Set lower and upper BC

for j $=1$:nx

$??$

end

$\%$ Set left and right BC

for i $=1$:nz

$??$

end

time $=0$;

for n$=1$:nt

$\%$ Compute rhs

rhs $=$ zeros$($nx$*$nz$,1)$;

for i $=1$:nz

for j $=1$:nx

ii $=$ Number$($i$,$j$)$;

$??$

end

end

$\%$ Compute solution vector

Tnew$\_$vector $=$ A$\backslash$rhs;

$\%$ Create $2$D matrix from vector

Tnew $=$ Tnew$\_$vector$($Number$)$;

T $=$ Tnew;

time $=$ time$+$dt;

$\%$ Plot solution every $50$ timesteps

if $($mod$($n$,10)==0)$

figure$(1),$ clf

pcolor$($x$2$d$/1$e$3,$z$2$d$/1$e$3,$Tnew$)$; shading interp, colorbar

hold on

contour$($x$2$d$/1$e$3,$z$2$d$/1$e$3,$Tnew,$[100$:$100$:$1500],\text{'}$k$\text{'})$;

xlabel$(\text{'}$x $[$km$]\text{'})$

ylabel$(\text{'}$z $[$km$]\text{'})$

zlabel$(\text{'}$Temperature $[^$oC$]\text{'})$

title$([\text{'}$Temperature evolution after $\text{'},$num$2$str$($time$/$year$/1$e$6),\text{'}$ Myrs$\text{'}])$

drawnow

end

end