diffusion equation in matlab

Hi, I have a pressure diffusion equation on a quadratic boundary. I have write the following code to solve it, the pressure should increase with time as we have injection in one side, and constant pressure other side. top and bottom side have isolated. but the code works only when length of medium is so small(<1). can anybody tell me how can I solve it for large length?

%B*DP/Dt+A*(DP/Dx+DP/Dy)+f(t)=0
%f(t)=-2.2/(t-0.1)^0.5
%General input data
clear

t0=0.1;
Sh=20;
G=20000;
nu=0.25;
n=0.716;
E=2*G*(1+nu);
K=1.2;
q_inj=0.067;
h=100;
h_tip=2500;
h_end=2600;
t_max=100;
N_t=60;
N_x=20;
N_y=30;
t = linspace(0.1,t_max,N_t);
y=linspace(-50,50,N_y);
h_fracture=linspace(h_tip,h_end,N_y);
B_y=((1-nu)/G)*(h^2-4*y.^2).^0.5;
B_y_mean=mean(B_y);
%Wellbore pressure calculation
vis=0.85;
den_f=1900;
den_p=2700;
prop_con=0.35;
den_eff=(1-prop_con)*den_f+prop_con*den_p;
D_t=1.28;
Re=4*q_inj*den_f/(pi*vis*D_t);
Re_w=((1+3*n)/4*n)*Re;
X_P1=log10(Re_w)/log10(exp(1));
X_P=log10(X_P1)/log10(exp(1));
f=exp(28.135+(-29.379+(8.2405-0.86227*X_P)*X_P)*X_P);
P_Frac=(-2*f*vis^2/(D_t*den_f)+den_f*9.81*0.001)*h_fracture*0.001;
w0_t=B_y'*(P_Frac-Sh);
w0=mean(mean(w0_t))
P_ave=mean(P_Frac)
B_y_mean=w0/P_ave
E_prin=E/(1-nu^2);
delt_t=(t_max-t0)/N_t;
Lf0=50;

delt_t=(t_max/(N_t));
delt_y=h/N_y;
C_Leak=(9.84*10^-6);
%pde geometry
A0=-w0^3/(12*vis);

rect=[3,4,0,Lf0,Lf0,0,0,0,100,100]';
gd=[rect];
ns=char('rect');
ns=ns';
sf='(rect)';
g=decsg(gd,sf,ns);
model=createpde;
geometryFromEdges(model,g);
generateMesh(model,'jiggle','on')
p = model.Mesh.Nodes;
pdegplot(model,'EdgeLabels','on')

%initial condition
u0 = Sh;
pdegplot(model,'EdgeLabels','on')
tlist = linspace(t0,t_max,N_t);
c=-A0*den_eff;
d=den_eff*B_y_mean;
f='-2.2./((t).^0.5)';
a=0;
applyBoundaryCondition(model,'Edge',4,'g',q_inj*den_eff/1000,'q',0)    
applyBoundaryCondition(model,'Edge',[1,3],'g',0,'q',0)
applyBoundaryCondition(model,'Edge',2,'u',Sh)
u1 = parabolic(u0,tlist,model,c,a,f,d);
sizeu1=size(u1);
nx=sizeu1(1);
x_w=linspace(0,Lf0,nx);
y_w=linspace(-h/2,h/2,nx);
B_y_w=((1-nu)/G)*(h^2-4*(p(2,:)-h/2).^2).^0.5;
w11=(u1(:,1)'-Sh).*B_y_w;
w12=(u1(:,N_t)'-Sh).*B_y_w;
max1=max(max(w12));
min1=min(min(w11));
figure1=figure
for ii=1:N_t
    hold off
   t=tlist(ii);

      figure1=pdeplot(model,'xydata',u1(:,ii));
      colormap default
      minu1=min(min(u1));
      maxu1=max(max(u1));
      caxis([minu1,maxu1])
  title(['t= ',num2str(t),'  sec'])
      pause(0.01)
      if ii<N_t
      delete(figure1)
      end
  end

figure2=figure
for ii=1:N_t
    hold off
    t=tlist(ii);
    w1=B_y_w'.*(u1(:,ii)-Sh);
    figure2=pdeplot(model,'xydata',w1);
    colormap default
    title(['t= ',num2str(t),'  sec'])
    caxis([min1,max1])
    colorbar
    shading interp
    pause(0.01)
    if ii<N_t
    delete(figure2)
    end
end

[Walter Roberson]

What happens if you use a larger medium?

Which variable corresponds to medium?