Help me,master-hands!

I am working for a week to solve this eqution by matlab,but I am lose.I don't know how to process the boundary condition.

who can tell me some code for solving the equation like this?
View attachment Microsoft Word (2).doc

Thanks a million !!!

 

yes, I know these are the boundary conditions.I mean that I don't know how to express them in the matlab program .

but thanks for your reply.

 

Thank you so much! you a so kind!

I wrote some code this afternoon.but i think there is some thing wrong

r is the longitudinal coordinates, t is the x-coordinate:

dr = 0.01; %each depth step
Nr = 100; % Choose the number of depth steps
Nt = 40; % Choose the number of time steps
dt = (10*365*24*60*60)/Nt; %Length of each time step in seconds

T = zeros(Nr+1,Nt+1); %Create temperature matrix with Nz+1 rows, and Nt+1 columns

time = [0:120/Nt:120]; %120 months in 10 years
T(:,1) = 1; %Set initial condition

maxiter = 500
for iter = 1:maxiter
%Tlast = T; %Save the last guess
%T(:,1) = Tlast(:,end); %Initialize the temp at t=0 to the last temp
for i=2:Nt+1
T(1,i)=(1+dt/dr)*T(1,i-1)-(dt/dr)*T(2,i-1);

for j=2:Nr

depth_2D = (1/(2*dr^2)*dr*(j-1)+(1-dr*(j-1))/(2*dr))*T(j-1,i-1)-(1/(dr^2)*dr*(j-1)+dr*(j-1))*T(j,i-1)+(1/(2*dr^2)*dr*(j-1)+(1-dr*(j-1))/(2*dr))*T(j+1,i-1);

T(j,i) = depth_2D*dt + T(j,i-1);

T(end,i) = T(end-1,i);
end
end
end

 

thank u so mush!but I think my difference method is implicit difference that is:
[tex] B_j^i = \frac{\sigma_j k}{h} B_{j-1}^{i+1} + \left(1-\frac{2 \sigma_j k}{h^2}-\frac{\mu_j k}{h}-r_j\right)B_j^{i+1} + \left(\frac{\sigma_j k}{h^2}+\frac{\mu_j k}{h}\right)B_{j+1}^{i+1},[/tex]

so:
[tex]B^i = \left( \begin{array}{cccccc} \alpha_2 & \beta_2 & \gamma_2 & 0 & 0 & \cdots \\ 0 & \alpha_3 & \beta_3 & \gamma_3 & 0 & \cdots \\ 0 & 0 & \alpha_4 & \beta_4 & \gamma_4 & \cdots \\ \vdots & & & & & \end{array}\right) B^{i+1} = \Lambda B^{i+1},[/tex]

maybe this is a equation system

 

because I use the implicit difference method.It is steady.So,maybe this is a equation system.

 

Thank you so much.
by the way,could I find the analyze solution of this equation by matlab?

 

thank you man! I am just back from norway.
I tried the code you give me,but it is not work.
The following is the code,could you help me to check it?

>>dr = 0.005;
Nr = 20;
dt = 0.1;
Nt = 20;
a=0.2339*0.0189;
b=0.2339;
delta=sqrt(0.0073);
u=zeros(Nr+1,Nt+1);
u(1,:)=1;

for j=3:Nr
alpha=(1/2)*delta^2*dr*(j-1);
mu=a-b*dr*(j-1);
rho=dt/(dr)^2;

A=diag((1-2*rho*alpha-(dt/dr)*mu-dt*dr*(j-1))*ones(Nr-2,1))-diag((rho*alpha+(dt/dr)*mu)*ones(Nr-3,1),1)-diag((rho*alpha)*ones(Nr-3,1),-1);
end
MatrixLambda = A; % Definition of Lambda.
for i=1:Nt
u(1,1) = u(1,1) - (a*dt/dr)*(u(2,1)-u(1,1)); % Left boundary (j = 0)
u(2,1) = rho*(1/2)*delta^2*dr*1* u(1,1) + 1-2*rho*(1/2)*delta^2*dr*1-(dt/dr)*(a-b*dr*1)-dt*dr*1* u(2,1) + rho*(1/2)*delta^2*dr*1+(dt/dr)*(a-b*dr*1)* u(3,1); % j = 1
u(3:Nr,1) = MatrixLambda*u(3:Nr,1); % j = 2, ..., n
u(Nr+1,1) = u(Nr,1); % Right boundary
V(:,i) = u(1:Nr+1,1); % Matrix V with the values of u^i in each column
end