Fdtxjr

$$u_t = 0.25(u_xx+u_yy) -8x(x^2+y^2), qquad 0 leq x,y leq 1, qquad 0 leq t leq 0.5,$$
$$u(0,y,t) = y^2 - y, qquad u(1,y,t) = y^2 + y + 1,$$
$$u(x,0,t) = x, qquad u(x,1,t) = 3x,$$
$$u(x,y,0) = y^" + 2xy + x -y.$$

I don't know how to implement the nonzero boundary conditions or transform the problem into a problem with zero boundary conditions.

Here is my attempt:

function [u,x,y,t]=q2(T,K,N,M)
% solves the 2D heat equation u_t=K u_xx+ f(x,y,t) in the unit square 
% and on the interval [0,T] with zero boundary condition by ADI method.
%K is the (heat conduction) coefficient in the heat equation.
%N is the number of subintervals in [0,1].
%M is the number of time steps.
h=1/N;
tau=T/M;
gamma=tau*K/h^2;
u=zeros(N+1,N+1,M+1);
u1=zeros(N+1,N+1);
alpha=zeros(N,1);
beta=zeros(N,1);
% grid points
x=(0:N)*h;
y=x;
t=(0:M)*tau;
%initial condition
for k=2:N
 for j=2:N
 u(k,j,1)=gg(x(k),y(j));
 end
end
% ADI method using the double-sweep method
for k=1:N-1
 alpha(k+1)=gamma/(2*(1+gamma)-gamma*alpha(k));
end
for n=2:M+1
 for j=2:N
 for k=2:N
 F=-0.5*gamma*(u(k,j-1,n-1)+u(k,j+1,n-1))-(1-gamma)*u(k,j,n-1)...
 -0.25*tau*(ff(x(k),y(j),t(n))+ff(x(k),y(j),t(n-1)));
 beta(k)=(beta(k-1)*gamma-2*F)/(2*(1+gamma)-gamma*alpha(k-1));
 end
 for k=N:(-1):2
 u1(k,j)=alpha(k)*u1(k+1,j)+beta(k);
 end
 end
 for k=2:N
 for j=2:N
 F=-0.5*gamma*(u1(k-1,j)+u1(k+1,j))-(1-gamma)*u1(k,j)...
 -0.25*tau*(ff(x(k),y(j),t(n))+ff(x(k),y(j),t(n-1)));
 beta(j)=(beta(j-1)*gamma-2*F)/(2*(1+gamma)-gamma*alpha(j-1));
 end
 for j=N:(-1):2
 u(k,j,n)=alpha(j)*u(k,j+1,n)+beta(j);
 end
 end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function f(x,y,t) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f=ff(x,y,t)
f=-8*x*(x^2+y^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function g(x,y) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function g=gg(x,y)
g=y^2+2*x*y+x-y;

For the general boundary conditions

beginalign
u(0,y,t) &= a_1(y) \
u(1,y,t) &= b_1(y) \
u(x,0,t) &= a_2(x) \
u(x,1,t) &= b_2(x)
endalign

You can write $u(x,y,t) = w_1(x,y) + w_2(x,y) + v(x,y,t)$ such that $v(x,y,t)$ is homogeneous and

beginarrayll
beginaligned
w_1(0,y) &= a_1(y) \
w_1(1,y) &= a_2(y)
endaligned &&&
beginaligned
w_2(x,0) &= a_2(x) - w_1(x,0) \
w_2(x,1) &= b_2(x) - w_1(x,1)
endaligned
endarray

The standard approach is to set
beginalign
w_1(x,y) &= (1-x)a_1(y) + xb_1(y) \
w_2(x,y) &= (1-y)[a_2(x) - (1-x)a_1(0) - xb_1(0)] + y[b_2(x) - (1-x)a_1(1) - xb_1(1)]
endalign

For this specific boundary problem, the simplified solution is

$$ w(x,y) = w_1(x,y) + w_2(x,y) = y^2 - y + x + 2xy $$

Note that this approach works only when the boundary is continuous at the "sharp points", i.e. $a_1(0)=a_2(0)$, $b_1(0)=a_2(1)$, $a_1(1)=b_2(0)$, $b_1(1)=b_2(1)$