Index exceed matrix dimensions.

% code for PTPB

Cv=0.5; %in m2/day

U0=100; %initial consolidation

U1s=0; %excess pore water pressure at top

U2s=0; %excess pore waler press at bottom

t=30; %in days

z=10; %total thickness in m

n=10; %no.of layers

dz=z/n; %length of interval

dt=1;%length of time interval

lambda= Cv*dt/dz^2;

m=t/dt;

U1=zeros(n,m);

U2=zeros(n,m);

U=zeros(n,m);

U1(2:n-1,1)=100;

U2(2:n-1,1)=100;

j=1;

for k=dt:dt:29*dt

j=j+1;

for i=2:n-1

U1(i,j)=(1-lambda)*U1(i,j-1)/(1+lambda)+(lambda)*U1(i+1,j-1)/(1+lambda)+(lambda)*U1(i-1,j)/(1+lambda);

end

end

for k=dt:dt:29*dt

j=j+1;

for i=n-1:-1:2

U2(i,j)=(1-lambda)*U2(i,j-1)/(1+lambda)+(lambda)*U2(i-1,j-1)/(1+lambda)+(lambda)*U2(i+1,j)/(1+lambda);

end

end

U=0.5*(U1+U2);

Answer 1

[Alan Stevens]

Do you just need

j = 1;

before the last

for k=dt:dt:29*dt

loop?

[Tarun Pandey]

can you help me how to plot this?

[Alan Stevens]

What exactly do you want to plot?

[Tarun Pandey]

U vs t (isochrones)

[Alan Stevens]

Do you mean something like this at the end of your code

[x, y] = meshgrid(0:dt:29*dt,1:n);
figure
surf(x,y,U1)
xlabel('x'),ylabel('y'),zlabel('U1')

and similarly for U2.

I should probably have used t instead of x!

[Alan Stevens]

Or more like this

figure
plot(0:dt:29*dt,U1),grid

[Tarun Pandey]

It worked thanku sir..