Index in position 2 exceeds array bounds (must not exceed 50).

1 view (last 30 days)
Dursman Mchabe
Dursman Mchabe on 19 Nov 2020
Commented: Dursman Mchabe on 19 Nov 2020
Hello everyone
I am trying to run the code below, however I get the error message:
logical
0
ans =
logical
0
Index in position 2 exceeds array bounds (must not exceed 50).
Error in Example (line 1348)
4/W .* A.*d2wdx2(i,j) - rr(i,j).*P(H.*I.*x.*L.*x - kb.*J.*x.*K.*x)==0; %Equation 36
>>
The code in question is
close all
clear all
clc
% Define number of rows and columns
Nt = 1000; % number of rows (for temporal steps)
Nx = 1000; % number of columns (for spatial steps)
% Define actual time and size
tMin = 1e-23; % minimum time
tMax = 100; % maximum time
xMin = 2; % minimum space
xMax = 100; % maximum space
% Define steps
t = linspace(tMin,tMax,Nt); % time steps
x = linspace(xMin,xMax,Nx); % spce steps
% Define actual stepping
delta_t = t(2) - t(1); % actual temporal stepping
delta_x = x(2) - x(1); % actual spacial spacing
%% GLOBAL
global A B C DD E F a G ...
H I J K L M N O P ...
Q R S T U V W X Y f Z ...
AA BB CC EE FF GG HH ...
II JJ KK LL MM NN OO ...
PP QQ RR SS r rprime g1 v a2 a6 alpha ...
mu CI kMT vfinal vinitial kb ...
amCa amMg Ap KHCO3 KH2O n TT UU VV WW XX YY ZZ AAA
log10(H) = 329.85 - 110.541*log10(T)-(17265.4./T); %(14)
%
M = (J.*rprime.*K.*rprime)/(I.*rprime.*L.*rprime); %(15)
%
KHCO3 = (J.*rprime.*U.*rprime)/(K.*rprime); %(16)
%
KH2O = (J.*rprime.*S.*rprime)/(L.*rprime); %(17)
%
rprime = 1 - (W/2*r); %(18)
%
% Z = - integral(AA.*n(v),vfinal,vinitial); %(21)
%
JJ = (kMT.*W)/B; %(23)
%
Ap = (pi* W^2)/4 == pi*((3*v)/(4*pi))^(2/3); %(24)
%
g1 = 10*exp(-7.84*v^2); %(26a)
%
a2 = -10*exp(-v); %(26b)
%
a6 = 25*tanh(10^3); %(26c)
%
alpha = -mu*v - ((6*g1^2 - (6*g1^2)/(18*v +1))/g1^2) + CI ; %(26d)
%
mu = 2.8;
%
CI = -4.0;
%
% n(v,t) = (1/3) + g1*exp(-(1/4)*lambertw((-a6*exp((4*sqrt(a2)*alpha)/(2*a2))))) + ...
% ((6*g1^2)*(exp(-(1/4)*lambertw((-a6*exp(4*sqrt(a2)*alpha)/(2*a2))+sqrt(a2)*alpha)^2)))/(18*v + 1); %(26)
% %
kb = H/M; %(39)
%
NN = amCa * U; %(43)
%
OO = amMg * U; %(44)
%
PP = amCa * amMg * (U)^2 ; %(45)
%-----------------------------------------------------------------------------------------------------------
%--------------------------------------------------------------------------
%% Pre-allocations of rows and columns
aa=zeros(Nt,Nx);
b=zeros(Nt,Nx);
c=zeros(Nt,Nx);
d=zeros(Nt,Nx);
e=zeros(Nt,Nx);
ff=zeros(Nt,Nx);
g=zeros(Nt,Nx);
h=zeros(Nt,Nx);
ii=zeros(Nt,Nx);
jj=zeros(Nt,Nx);
k=zeros(Nt,Nx);
ll=zeros(Nt,Nx);
m=zeros(Nt,Nx);
nn=zeros(Nt,Nx);
o=zeros(Nt,Nx);
p=zeros(Nt,Nx);
q=zeros(Nt,Nx);
rr=zeros(Nt,Nx);
s=zeros(Nt,Nx);
tt=zeros(Nt,Nx);
u=zeros(Nt,Nx);
v=zeros(Nt,Nx);
w=zeros(Nt,Nx);
xx=zeros(Nt,Nx);
yy=zeros(Nt,Nx);
zz=zeros(Nt,Nx);
aaa=zeros(Nt,Nx);
bb=zeros(Nt,Nx);
cc=zeros(Nt,Nx);
dd=zeros(Nt,Nx);
ee=zeros(Nt,Nx);
% -------------------------------------------------------------------------
A = 1;
B = 1;
C = 1;
DD = 1;
E = 1;
a = 1;
G = 1;
H = 1;
I = 1;
J = 1;
K = 1;
L = 1;
M = 1;
N = 1;
O = 1;
P = 1;
Q = 1;
R = 1;
S = 1;
T = 1;
U = 1;
V = 1;
W = 1;
X = a;
a = 1;
Y = f;
f = 1;
Z = 1;
AA = 1;
BB = 1;
CC = 1;
L = 1;
EE = 1;
FF = 0;
GG = 5e-6;
HH = 1500e-6;
II = 1;
JJ = 1;
KK = 1;
LL = 1;
MM = 1;
NN = 1;
OO = 1;
PP = 1;
QQ = 1;
RR = 1;
SS = 1;
TT =1;
UU = 1;
VV = 1;
WW = 1;
XX = 1;
YY = 1;
ZZ = 1;
AAA =1;
% Initial conditions (at t = tMin)
aa(1,:)=0;
b(1,:)=0;
c(1,:)=0;
d(1,:)=0;
e(1,:)=0;
ff(1,:)=0;
g(1,:)=0;
h(1,:)=0;
ii(1,:)=0;
jj(1,:)=0;
k(1,:)=0;
ll(1,:)=0;
m(1,:)=0;
nn(1,:)=0;
o(1,:)=0;
p(1,:)=0;
q(1,:)=0;
rr(1,:)=0;
s(1,:)=0;
tt(1,:)=0;
u(1,:)=0;
v(1,:)=0;
w(1,:)=0;
xx(1,:)=0;
yy(1,:)=0;
zz(1,:)=0;
aaa(1,:)=0;
bb(1,:)=0;
cc(1,:)=0;
dd(1,:)=0;
ee(1,:)=0;
% Second row to final row (time iterations)
for i = 2:Nt
% First column to final colun (space iterations)
for j = 2:Nx
% First column
if j==2
% DEFINING TERMS
% previous time point (i-1) and forward spatial point (j+1)
aa_im1_jp1=aa(i-1,j+1);
b_im1_jp1=b(i-1,j+1);
c_im1_jp1=c(i-1,j+1);
d_im1_jp1=d(i-1,j+1);
e_im1_jp1=e(i-1,j+1);
ff_im1_jp1=ff(i-1,j+1);
g_im1_jp1=g(i-1,j+1);
h_im1_jp1=h(i-1,j+1);
ii_im1_jp1=ii(i-1,j+1);
jj_im1_jp1=jj(i-1,j+1);
k_im1_jp1=k(i-1,j+1);
ll_im1_jp1=ll(i-1,j+1);
m_im1_jp1=m(i-1,j+1);
nn_im1_jp1=nn(i-1,j+1);
o_im1_jp1=o(i-1,j+1);
p_im1_jp1=p(i-1,j+1);
q_im1_jp1=q(i-1,j+1);
rr_im1_jp1=rr(i-1,j+1);
s_im1_jp1=s(i-1,j+1);
tt_im1_jp1=tt(i-1,j+1);
u_im1_jp1=u(i-1,j+1);
v_im1_jp1=v(i-1,j+1);
w_im1_jp1=w(i-1,j+1);
xx_im1_jp1=x(i-1,j+1);
yy_im1_jp1=yy(i-1,j+1);
zz_im1_jp1=zz(i-1,j+1);
aaa_im1_jp1=aaa(i-1,j+1);
bb_im1_jp1=bb(i-1,j+1);
cc_im1_jp1=cc(i-1,j+1);
dd_im1_jp1=dd(i-1,j+1);
ee_im1_jp1=ee(i-1,j+1);
% previous time point (i-1) and current spatial point (j)
aa_im1_j=aa(i-1,j);
b_im1_j=b(i-1,j);
c_im1_j=c(i-1,j);
d_im1_j=d(i-1,j);
e_im1_j=e(i-1,j);
ff_im1_j=ff(i-1,j);
g_im1_j=g(i-1,j);
h_im1_j=h(i-1,j);
ii_im1_j=ii(i-1,j);
jj_im1_j=jj(i-1,j);
k_im1_j=k(i-1,j);
ll_im1_j=ll(i-1,j);
m_im1_j=m(i-1,j);
nn_im1_j=nn(i-1,j);
o_im1_j=o(i-1,j);
p_im1_j=p(i-1,j);
q_im1_j=q(i-1,j);
rr_im1_j=rr(i-1,j);
s_im1_j=s(i-1,j);
tt_im1_j=tt(i-1,j);
u_im1_j=u(i-1,j);
v_im1_j=v(i-1,j);
w_im1_j=w(i-1,j);
xx_im1_j=x(i-1,j);
yy_im1_j=yy(i-1,j);
zz_im1_j=zz(i-1,j);
aaa_im1_j=aaa(i-1,j);
bb_im1_j=bb(i-1,j);
cc_im1_j=cc(i-1,j);
dd_im1_j=dd(i-1,j);
ee_im1_j=ee(i-1,j);
% previous time point (i-1) and previous spatial point (j-1)
aa_im1_jm1=aa(i-1,j-1);
b_im1_jm1=b(i-1,j-1);
c_im1_jm1=c(i-1,j-1);
d_im1_jm1=d(i-1,j-1);
e_im1_jm1=e(i-1,j-1);
ff_im1_jm1=ff(i-1,j-1);
g_im1_jm1=g(i-1,j-1);
h_im1_jm1=h(i-1,j-1);
ii_im1_jm1=ii(i-1,j-1);
jj_im1_jm1=jj(i-1,j-1);
k_im1_jm1=k(i-1,j-1);
ll_im1_jm1=ll(i-1,j-1);
m_im1_jm1=m(i-1,j-1);
nn_im1_jm1=nn(i-1,j-1);
o_im1_jm1=o(i-1,j-1);
p_im1_jm1=p(i-1,j-1);
q_im1_jm1=q(i-1,j-1);
rr_im1_jm1=rr(i-1,j-1);
s_im1_jm1=s(i-1,j-1);
tt_im1_jm1=tt(i-1,j-1);
u_im1_jm1=u(i-1,j-1);
v_im1_jm1=v(i-1,j-1);
w_im1_jm1=w(i-1,j-1);
xx_im1_jm1=x(i-1,j-1);
yy_im1_jm1=yy(i-1,j-1);
zz_im1_jm1=zz(i-1,j-1);
aaa_im1_jm1=aaa(i-1,j-1);
bb_im1_jm1=bb(i-1,j-1);
cc_im1_jm1=cc(i-1,j-1);
dd_im1_jm1=dd(i-1,j-1);
ee_im1_jm1=ee(i-1,j-1);
% current time point (i) and next spatial point (j+1)
aa_i_jp1=aa(i,j+1);
b_i_jp1=b(i,j+1);
c_i_jp1=c(i,j+1);
d_i_jp1=d(i,j+1);
e_i_jp1=e(i,j+1);
ff_i_jp1=ff(i,j+1);
g_i_jp1=g(i,j+1);
h_i_jp1=h(i,j+1);
ii_i_jp1=ii(i,j+1);
jj_i_jp1=jj(i,j+1);
k_i_jp1=k(i,j+1);
ll_i_jp1=ll(i,j+1);
m_i_jp1=m(i,j+1);
nn_i_jp1=nn(i,j+1);
o_i_jp1=o(i,j+1);
p_i_jp1=p(i,j+1);
q_i_jp1=q(i,j+1);
rr_i_jp1=rr(i,j+1);
s_i_jp1=s(i,j+1);
tt_i_jp1=tt(i,j+1);
u_i_jp1=u(i,j+1);
v_i_jp1=v(i,j+1);
w_i_jp1=w(i,j+1);
xx_i_jp1=xx(i,j+1);
yy_i_jp1=yy(i,j+1);
zz_i_jp1=zz(i,j+1);
aaa_i_jp1=aaa(i,j+1);
bb_i_jp1=bb(i,j+1);
cc_i_jp1=cc(i,j+1);
dd_i_jp1=dd(i,j+1);
ee_i_jp1=ee(i,j+1);
% current time point (i) and previous spatial point (j-1)
aa_i_jm1=aa(i,j-1);
b_i_jm1=b(i,j-1);
c_i_jm1=c(i,j-1);
d_i_jm1=d(i,j-1);
e_i_jm1=e(i,j-1);
ff_i_jm1=ff(i,j-1);
g_i_jm1=g(i,j-1);
h_i_jm1=h(i,j-1);
ii_i_jm1=ii(i,j-1);
jj_i_jm1=jj(i,j-1);
k_i_jm1=k(i,j-1);
ll_i_jm1=ll(i,j-1);
m_i_jm1=m(i,j-1);
nn_i_jm1=nn(i,j-1);
o_i_jm1=o(i,j-1);
p_i_jm1=p(i,j-1);
q_i_jm1=q(i,j-1);
rr_i_jm1=rr(i,j-1);
s_i_jm1=s(i,j-1);
tt_i_jm1=tt(i,j-1);
u_i_jm1=u(i,j-1);
v_i_jm1=v(i,j-1);
w_i_jm1=w(i,j-1);
xx_i_jm1=xx(i,j-1);
yy_i_jm1=yy(i,j-1);
zz_i_jm1=zz(i,j-1);
aaa_i_jm1=aaa(i,j-1);
bb_i_jm1=bb(i,j-1);
cc_i_jm1=cc(i,j-1);
dd_i_jm1=dd(i,j-1);
ee_i_jm1=ee(i,j-1);
% Discretization of space first order PDE terms using central difference method of lines
daadx(i,j)=(aa_i_jp1-aa_i_jm1)/(2*delta_x);
dbdx(i,j)=(b_i_jp1-b_i_jm1)/(2*delta_x);
dcdx(i,j)=(c_i_jp1-c_i_jm1)/(2*delta_x);
dddx(i,j)=(dd_i_jp1-dd_i_jm1)/(2*delta_x);
dedx(i,j)=(e_i_jp1-e_i_jm1)/(2*delta_x);
dffdx(i,j)=(ff_i_jp1-ff_i_jm1)/(2*delta_x);
dgdx(i,j)=(g_i_jp1-g_i_jm1)/(2*delta_x);
dhdx(i,j)=(h_i_jp1-h_i_jm1)/(2*delta_x);
diidx(i,j)=(ii_i_jp1-ii_i_jm1)/(2*delta_x);
djjdx(i,j)=(jj_i_jp1-jj_i_jm1)/(2*delta_x);
dkdx(i,j)=(k_i_jp1-k_i_jm1)/(2*delta_x);
dlldx(i,j)=(ll_i_jp1-ll_i_jm1)/(2*delta_x);
dmdx(i,j)=(m_i_jp1-m_i_jm1)/(2*delta_x);
dnndx(i,j)=(nn_i_jp1-nn_i_jm1)/(2*delta_x);
dodx(i,j)=(o_i_jp1-o_i_jm1)/(2*delta_x);
Wdx(i,j)=(p_i_jp1-p_i_jm1)/(2*delta_x);
dqdx(i,j)=(q_i_jp1-q_i_jm1)/(2*delta_x);
drrdx(i,j)=(rr_i_jp1-rr_i_jm1)/(2*delta_x);
dsdx(i,j)=(s_i_jp1-s_i_jm1)/(2*delta_x);
dttdx(i,j)=(tt_i_jp1-tt_i_jm1)/(2*delta_x);
dudx(i,j)=(u_i_jp1-u_i_jm1)/(2*delta_x);
dvdx(i,j)=(v_i_jp1-v_i_jm1)/(2*delta_x);
dwdx(i,j)=(w_i_jp1-w_i_jm1)/(2*delta_x);
dxxdx(i,j)=(xx_i_jp1-xx_i_jm1)/(2*delta_x);
dyydx(i,j)=(yy_i_jp1-yy_i_jm1)/(2*delta_x);
dzzdx(i,j)=(zz_i_jp1-zz_i_jm1)/(2*delta_x);
daaadx(i,j)=(aaa_i_jp1-aaa_i_jm1)/(2*delta_x);
dbbdx(i,j)=(bb_i_jp1-bb_i_jm1)/(2*delta_x);
dccdx(i,j)=(cc_i_jp1-cc_i_jm1)/(2*delta_x);
ddddx(i,j)=(dd_i_jp1-dd_i_jm1)/(2*delta_x);
deedx(i,j)=(ee_i_jp1-ee_i_jm1)/(2*delta_x);
% Discretization of space second order PDE terms using central difference method of lines
d2aadx2(i,j)=(aa_i_jp1-2.*aa(i,j)+aa_i_jm1)/(delta_x.^2);
d2bdx2(i,j)=(b_i_jp1-2.*b(i,j)+b_i_jm1)/(delta_x.^2);
d2cdx2(i,j)=(c_i_jp1-2.*c(i,j)+c_i_jm1)/(delta_x.^2);
d2dddx2(i,j)=(dd_i_jp1-2.*dd(i,j)+dd_i_jm1)/(delta_x.^2);
d2edx2(i,j)=(e_i_jp1-2.*e(i,j)+e_i_jm1)/(delta_x.^2);
d2ffdx2(i,j)=(ff_i_jp1-2.*ff(i,j)+ff_i_jm1)/(delta_x.^2);
d2gdx2(i,j)=(g_i_jp1-2.*g(i,j)+g_i_jm1)/(delta_x.^2);
d2hdx2(i,j)=(h_i_jp1-2.*h(i,j)+h_i_jm1)/(delta_x.^2);
d2iidx2(i,j)=(ii_i_jp1-2.*ii(i,j)+ii_i_jm1)/(delta_x.^2);
d2jjdx2(i,j)=(jj_i_jp1-2.*jj(i,j)+jj_i_jm1)/(delta_x.^2);
d2kdx2(i,j)=(k_i_jp1-2.*k(i,j)+k_i_jm1)/(delta_x.^2);
d2lldx2(i,j)=(ll_i_jp1-2.*ll(i,j)+ll_i_jm1)/(delta_x.^2);
d2mdx2(i,j)=(m_i_jp1-2.*m(i,j)+m_i_jm1)/(delta_x.^2);
d2nndx2(i,j)=(nn_i_jp1-2.*nn(i,j)+nn_i_jm1)/(delta_x.^2);
d2odx2(i,j)=(o_i_jp1-2.*o(i,j)+o_i_jm1)/(delta_x.^2);
d2pdx2(i,j)=(p_i_jp1-2.*p(i,j)+p_i_jm1)/(delta_x.^2);
d2qdx2(i,j)=(q_i_jp1-2.*q(i,j)+q_i_jm1)/(delta_x.^2);
d2rrdx2(i,j)=(rr_i_jp1-2.*rr(i,j)+rr_i_jm1)/(delta_x.^2);
d2sdx2(i,j)=(s_i_jp1-2.*s(i,j)+s_i_jm1)/(delta_x.^2);
d2ttdx2(i,j)=(tt_i_jp1-2.*tt(i,j)+tt_i_jm1)/(delta_x.^2);
d2udx2(i,j)=(u_i_jp1-2.*u(i,j)+u_i_jm1)/(delta_x.^2);
d2vdx2(i,j)=(v_i_jp1-2.*v(i,j)+v_i_jm1)/(delta_x.^2);
d2wdx2(i,j)=(w_i_jp1-2.*w(i,j)+w_i_jm1)/(delta_x.^2);
d2xxdx2(i,j)=(xx_i_jp1-2.*xx(i,j)+xx_i_jm1)/(delta_x.^2);
d2yydx2(i,j)=(yy_i_jp1-2.*yy(i,j)+yy_i_jm1)/(delta_x.^2);
d2zzdx2(i,j)=(zz_i_jp1-2.*zz(i,j)+zz_i_jm1)/(delta_x.^2);
d2aaadx2(i,j)=(aaa_i_jp1-2.*aaa(i,j)+aaa_i_jm1)/(delta_x.^2);
d2bbdx2(i,j)=(bb_i_jp1-2.*bb(i,j)+bb_i_jm1)/(delta_x.^2);
d2ccdx2(i,j)=(cc_i_jp1-2.*cc(i,j)+cc_i_jm1)/(delta_x.^2);
d2dddx2(i,j)=(dd_i_jp1-2.*dd(i,j)+dd_i_jm1)/(delta_x.^2);
d2eedx2(i,j)=(ee_i_jp1-2.*ee(i,j)+ee_i_jm1)/(delta_x.^2);
% Substituting equations as they appear in the paper
A.*QQ.*daadx(i,j)+ G.*a.*(b(i,j)-c(i,j))==0; % Equation
B.*dddx(i,j)==0; % Equation 6
C.*dedx(i,j)==0; % Equation 7
E.*dffdx(i,j) + TT.*dgdx(i,j)==0; % Equation 8
UU.*B.*dddx(i,j)+ VV.*C.*dedx(i,j) + WW.*AAA.*diidx(i,j)...
- XX.*DD.*djjdx(i,j) - YY.*E.*dffdx(i,j) - ZZ.* TT.*dgdx(i,j) == 0; % Equation 9
%==================================================================================================================================================================
elseif j>=1 & j<=50
% previous time point (i-1) and forward spatial point (j+1)
aa_im1_jp1=aa(i-1,j+1);
b_im1_jp1=b(i-1,j+1);
c_im1_jp1=c(i-1,j+1);
d_im1_jp1=d(i-1,j+1);
e_im1_jp1=e(i-1,j+1);
ff_im1_jp1=ff(i-1,j+1);
g_im1_jp1=g(i-1,j+1);
h_im1_jp1=h(i-1,j+1);
ii_im1_jp1=ii(i-1,j+1);
jj_im1_jp1=jj(i-1,j+1);
k_im1_jp1=k(i-1,j+1);
ll_im1_jp1=ll(i-1,j+1);
m_im1_jp1=m(i-1,j+1);
nn_im1_jp1=nn(i-1,j+1);
o_im1_jp1=o(i-1,j+1);
p_im1_jp1=p(i-1,j+1);
q_im1_jp1=q(i-1,j+1);
rr_im1_jp1=rr(i-1,j+1);
s_im1_jp1=s(i-1,j+1);
tt_im1_jp1=tt(i-1,j+1);
u_im1_jp1=u(i-1,j+1);
v_im1_jp1=v(i-1,j+1);
w_im1_jp1=w(i-1,j+1);
xx_im1_jp1=x(i-1,j+1);
yy_im1_jp1=yy(i-1,j+1);
zz_im1_jp1=zz(i-1,j+1);
aaa_im1_jp1=aaa(i-1,j+1);
bb_im1_jp1=bb(i-1,j+1);
cc_im1_jp1=cc(i-1,j+1);
dd_im1_jp1=dd(i-1,j+1);
ee_im1_jp1=ee(i-1,j+1);
% previous time point (i-1) and current spatial point (j)
aa_im1_j=aa(i-1,j);
b_im1_j=b(i-1,j);
c_im1_j=c(i-1,j);
d_im1_j=d(i-1,j);
e_im1_j=e(i-1,j);
ff_im1_j=ff(i-1,j);
g_im1_j=g(i-1,j);
h_im1_j=h(i-1,j);
ii_im1_j=ii(i-1,j);
jj_im1_j=jj(i-1,j);
k_im1_j=k(i-1,j);
ll_im1_j=ll(i-1,j);
m_im1_j=m(i-1,j);
nn_im1_j=nn(i-1,j);
o_im1_j=o(i-1,j);
p_im1_j=p(i-1,j);
q_im1_j=q(i-1,j);
rr_im1_j=rr(i-1,j);
s_im1_j=s(i-1,j);
tt_im1_j=tt(i-1,j);
u_im1_j=u(i-1,j);
v_im1_j=v(i-1,j);
w_im1_j=w(i-1,j);
xx_im1_j=x(i-1,j);
yy_im1_j=yy(i-1,j);
zz_im1_j=zz(i-1,j);
aaa_im1_j=aaa(i-1,j);
bb_im1_j=bb(i-1,j);
cc_im1_j=cc(i-1,j);
dd_im1_j=dd(i-1,j);
ee_im1_j=ee(i-1,j);
% previous time point (i-1) and previous spatial point (j-1)
aa_im1_jm1=aa(i-1,j-1);
b_im1_jm1=b(i-1,j-1);
c_im1_jm1=c(i-1,j-1);
d_im1_jm1=d(i-1,j-1);
e_im1_jm1=e(i-1,j-1);
ff_im1_jm1=ff(i-1,j-1);
g_im1_jm1=g(i-1,j-1);
h_im1_jm1=h(i-1,j-1);
ii_im1_jm1=ii(i-1,j-1);
jj_im1_jm1=jj(i-1,j-1);
k_im1_jm1=k(i-1,j-1);
ll_im1_jm1=ll(i-1,j-1);
m_im1_jm1=m(i-1,j-1);
nn_im1_jm1=nn(i-1,j-1);
o_im1_jm1=o(i-1,j-1);
p_im1_jm1=p(i-1,j-1);
q_im1_jm1=q(i-1,j-1);
rr_im1_jm1=rr(i-1,j-1);
s_im1_jm1=s(i-1,j-1);
tt_im1_jm1=tt(i-1,j-1);
u_im1_jm1=u(i-1,j-1);
v_im1_jm1=v(i-1,j-1);
w_im1_jm1=w(i-1,j-1);
xx_im1_jm1=x(i-1,j-1);
yy_im1_jm1=yy(i-1,j-1);
zz_im1_jm1=zz(i-1,j-1);
aaa_im1_jm1=aaa(i-1,j-1);
bb_im1_jm1=bb(i-1,j-1);
cc_im1_jm1=cc(i-1,j-1);
dd_im1_jm1=dd(i-1,j-1);
ee_im1_jm1=ee(i-1,j-1);
% current time point (i) and next spatial point (j+1)
aa_i_jp1=aa(i,j+1);
b_i_jp1=b(i,j+1);
c_i_jp1=c(i,j+1);
d_i_jp1=d(i,j+1);
e_i_jp1=e(i,j+1);
ff_i_jp1=ff(i,j+1);
g_i_jp1=g(i,j+1);
h_i_jp1=h(i,j+1);
ii_i_jp1=ii(i,j+1);
jj_i_jp1=jj(i,j+1);
k_i_jp1=k(i,j+1);
ll_i_jp1=ll(i,j+1);
m_i_jp1=m(i,j+1);
nn_i_jp1=nn(i,j+1);
o_i_jp1=o(i,j+1);
p_i_jp1=p(i,j+1);
q_i_jp1=q(i,j+1);
rr_i_jp1=rr(i,j+1);
s_i_jp1=s(i,j+1);
tt_i_jp1=tt(i,j+1);
u_i_jp1=u(i,j+1);
v_i_jp1=v(i,j+1);
w_i_jp1=w(i,j+1);
xx_i_jp1=xx(i,j+1);
yy_i_jp1=yy(i,j+1);
zz_i_jp1=zz(i,j+1);
aaa_i_jp1=aaa(i,j+1);
bb_i_jp1=bb(i,j+1);
cc_i_jp1=cc(i,j+1);
dd_i_jp1=dd(i,j+1);
ee_i_jp1=ee(i,j+1);
% current time point (i) and previous spatial point (j-1)
aa_i_jm1=aa(i,j-1);
b_i_jm1=b(i,j-1);
c_i_jm1=c(i,j-1);
d_i_jm1=d(i,j-1);
e_i_jm1=e(i,j-1);
ff_i_jm1=ff(i,j-1);
g_i_jm1=g(i,j-1);
h_i_jm1=h(i,j-1);
ii_i_jm1=ii(i,j-1);
jj_i_jm1=jj(i,j-1);
k_i_jm1=k(i,j-1);
ll_i_jm1=ll(i,j-1);
m_i_jm1=m(i,j-1);
nn_i_jm1=nn(i,j-1);
o_i_jm1=o(i,j-1);
p_i_jm1=p(i,j-1);
q_i_jm1=q(i,j-1);
rr_i_jm1=rr(i,j-1);
s_i_jm1=s(i,j-1);
tt_i_jm1=tt(i,j-1);
u_i_jm1=u(i,j-1);
v_i_jm1=v(i,j-1);
w_i_jm1=w(i,j-1);
xx_i_jm1=xx(i,j-1);
yy_i_jm1=yy(i,j-1);
zz_i_jm1=zz(i,j-1);
aaa_i_jm1=aaa(i,j-1);
bb_i_jm1=bb(i,j-1);
cc_i_jm1=cc(i,j-1);
dd_i_jm1=dd(i,j-1);
ee_i_jm1=ee(i,j-1);
% Discretization of space first order PDE terms using central difference method of lines
daadx(i,j)=(aa_i_jp1-aa_i_jm1)/(2*delta_x);
dbdx(i,j)=(b_i_jp1-b_i_jm1)/(2*delta_x);
dcdx(i,j)=(c_i_jp1-c_i_jm1)/(2*delta_x);
dddx(i,j)=(dd_i_jp1-dd_i_jm1)/(2*delta_x);
dedx(i,j)=(e_i_jp1-e_i_jm1)/(2*delta_x);
dffdx(i,j)=(ff_i_jp1-ff_i_jm1)/(2*delta_x);
dgdx(i,j)=(g_i_jp1-g_i_jm1)/(2*delta_x);
dhdx(i,j)=(h_i_jp1-h_i_jm1)/(2*delta_x);
diidx(i,j)=(ii_i_jp1-ii_i_jm1)/(2*delta_x);
djjdx(i,j)=(jj_i_jp1-jj_i_jm1)/(2*delta_x);
dkdx(i,j)=(k_i_jp1-k_i_jm1)/(2*delta_x);
dlldx(i,j)=(ll_i_jp1-ll_i_jm1)/(2*delta_x);
dmdx(i,j)=(m_i_jp1-m_i_jm1)/(2*delta_x);
dnndx(i,j)=(nn_i_jp1-nn_i_jm1)/(2*delta_x);
dodx(i,j)=(o_i_jp1-o_i_jm1)/(2*delta_x);
Wdx(i,j)=(p_i_jp1-p_i_jm1)/(2*delta_x);
dqdx(i,j)=(q_i_jp1-q_i_jm1)/(2*delta_x);
drrdx(i,j)=(rr_i_jp1-rr_i_jm1)/(2*delta_x);
dsdx(i,j)=(s_i_jp1-s_i_jm1)/(2*delta_x);
dttdx(i,j)=(tt_i_jp1-tt_i_jm1)/(2*delta_x);
dudx(i,j)=(u_i_jp1-u_i_jm1)/(2*delta_x);
dvdx(i,j)=(v_i_jp1-v_i_jm1)/(2*delta_x);
dwdx(i,j)=(w_i_jp1-w_i_jm1)/(2*delta_x);
dxxdx(i,j)=(xx_i_jp1-xx_i_jm1)/(2*delta_x);
dyydx(i,j)=(yy_i_jp1-yy_i_jm1)/(2*delta_x);
dzzdx(i,j)=(zz_i_jp1-zz_i_jm1)/(2*delta_x);
daaadx(i,j)=(aaa_i_jp1-aaa_i_jm1)/(2*delta_x);
dbbdx(i,j)=(bb_i_jp1-bb_i_jm1)/(2*delta_x);
dccdx(i,j)=(cc_i_jp1-cc_i_jm1)/(2*delta_x);
ddddx(i,j)=(dd_i_jp1-dd_i_jm1)/(2*delta_x);
deedx(i,j)=(ee_i_jp1-ee_i_jm1)/(2*delta_x);
% Discretization of space second order PDE terms using central difference method of lines
d2aadx2(i,j)=(aa_i_jp1-2.*aa(i,j)+aa_i_jm1)/(delta_x.^2);
d2bdx2(i,j)=(b_i_jp1-2.*b(i,j)+b_i_jm1)/(delta_x.^2);
d2cdx2(i,j)=(c_i_jp1-2.*c(i,j)+c_i_jm1)/(delta_x.^2);
d2dddx2(i,j)=(dd_i_jp1-2.*dd(i,j)+dd_i_jm1)/(delta_x.^2);
d2edx2(i,j)=(e_i_jp1-2.*e(i,j)+e_i_jm1)/(delta_x.^2);
d2ffdx2(i,j)=(ff_i_jp1-2.*ff(i,j)+ff_i_jm1)/(delta_x.^2);
d2gdx2(i,j)=(g_i_jp1-2.*g(i,j)+g_i_jm1)/(delta_x.^2);
d2hdx2(i,j)=(h_i_jp1-2.*h(i,j)+h_i_jm1)/(delta_x.^2);
d2iidx2(i,j)=(ii_i_jp1-2.*ii(i,j)+ii_i_jm1)/(delta_x.^2);
d2jjdx2(i,j)=(jj_i_jp1-2.*jj(i,j)+jj_i_jm1)/(delta_x.^2);
d2kdx2(i,j)=(k_i_jp1-2.*k(i,j)+k_i_jm1)/(delta_x.^2);
d2lldx2(i,j)=(ll_i_jp1-2.*ll(i,j)+ll_i_jm1)/(delta_x.^2);
d2mdx2(i,j)=(m_i_jp1-2.*m(i,j)+m_i_jm1)/(delta_x.^2);
d2nndx2(i,j)=(nn_i_jp1-2.*nn(i,j)+nn_i_jm1)/(delta_x.^2);
d2odx2(i,j)=(o_i_jp1-2.*o(i,j)+o_i_jm1)/(delta_x.^2);
d2pdx2(i,j)=(p_i_jp1-2.*p(i,j)+p_i_jm1)/(delta_x.^2);
d2qdx2(i,j)=(q_i_jp1-2.*q(i,j)+q_i_jm1)/(delta_x.^2);
d2rrdx2(i,j)=(rr_i_jp1-2.*rr(i,j)+rr_i_jm1)/(delta_x.^2);
d2sdx2(i,j)=(s_i_jp1-2.*s(i,j)+s_i_jm1)/(delta_x.^2);
d2ttdx2(i,j)=(tt_i_jp1-2.*tt(i,j)+tt_i_jm1)/(delta_x.^2);
d2udx2(i,j)=(u_i_jp1-2.*u(i,j)+u_i_jm1)/(delta_x.^2);
d2vdx2(i,j)=(v_i_jp1-2.*v(i,j)+v_i_jm1)/(delta_x.^2);
d2wdx2(i,j)=(w_i_jp1-2.*w(i,j)+w_i_jm1)/(delta_x.^2);
d2xxdx2(i,j)=(xx_i_jp1-2.*xx(i,j)+xx_i_jm1)/(delta_x.^2);
d2yydx2(i,j)=(yy_i_jp1-2.*yy(i,j)+yy_i_jm1)/(delta_x.^2);
d2zzdx2(i,j)=(zz_i_jp1-2.*zz(i,j)+zz_i_jm1)/(delta_x.^2);
d2aaadx2(i,j)=(aaa_i_jp1-2.*aaa(i,j)+aaa_i_jm1)/(delta_x.^2);
d2bbdx2(i,j)=(bb_i_jp1-2.*bb(i,j)+bb_i_jm1)/(delta_x.^2);
d2ccdx2(i,j)=(cc_i_jp1-2.*cc(i,j)+cc_i_jm1)/(delta_x.^2);
d2dddx2(i,j)=(dd_i_jp1-2.*dd(i,j)+dd_i_jm1)/(delta_x.^2);
d2eedx2(i,j)=(ee_i_jp1-2.*ee(i,j)+ee_i_jm1)/(delta_x.^2);
% Discretization of time first order PDE terms using central difference method of lines
daadt(i,j)=(aa_i_jp1-aa_i_jm1)/(2*delta_t);
dbdt(i,j)=(b_i_jp1-b_i_jm1)/(2*delta_t);
dcdt(i,j)=(c_i_jp1-c_i_jm1)/(2*delta_t);
dddt(i,j)=(dd_i_jp1-dd_i_jm1)/(2*delta_t);
dedt(i,j)=(e_i_jp1-e_i_jm1)/(2*delta_t);
dffdt(i,j)=(ff_i_jp1-ff_i_jm1)/(2*delta_t);
dgdt(i,j)=(g_i_jp1-g_i_jm1)/(2*delta_t);
dhdt(i,j)=(h_i_jp1-h_i_jm1)/(2*delta_t);
diidt(i,j)=(ii_i_jp1-ii_i_jm1)/(2*delta_t);
djjdt(i,j)=(jj_i_jp1-jj_i_jm1)/(2*delta_t);
dkdt(i,j)=(k_i_jp1-k_i_jm1)/(2*delta_t);
dlldt(i,j)=(ll_i_jp1-ll_i_jm1)/(2*delta_t);
dmdt(i,j)=(m_i_jp1-m_i_jm1)/(2*delta_t);
dnndt(i,j)=(nn_i_jp1-nn_i_jm1)/(2*delta_t);
dodt(i,j)=(o_i_jp1-o_i_jm1)/(2*delta_t);
Wdt(i,j)=(p_i_jp1-p_i_jm1)/(2*delta_t);
dqdt(i,j)=(q_i_jp1-q_i_jm1)/(2*delta_t);
drrdt(i,j)=(rr_i_jp1-rr_i_jm1)/(2*delta_t);
dsdt(i,j)=(s_i_jp1-s_i_jm1)/(2*delta_t);
dttdt(i,j)=(tt_i_jp1-tt_i_jm1)/(2*delta_t);
dudt(i,j)=(u_i_jp1-u_i_jm1)/(2*delta_t);
dvdt(i,j)=(v_i_jp1-v_i_jm1)/(2*delta_t);
dwdt(i,j)=(w_i_jp1-w_i_jm1)/(2*delta_t);
dxxdt(i,j)=(xx_i_jp1-xx_i_jm1)/(2*delta_t);
dyydt(i,j)=(yy_i_jp1-yy_i_jm1)/(2*delta_t);
dzzdt(i,j)=(zz_i_jp1-zz_i_jm1)/(2*delta_t);
daaadt(i,j)=(aaa_i_jp1-aaa_i_jm1)/(2*delta_t);
dbbdt(i,j)=(bb_i_jp1-bb_i_jm1)/(2*delta_t);
dccdt(i,j)=(cc_i_jp1-cc_i_jm1)/(2*delta_t);
ddddt(i,j)=(dd_i_jp1-dd_i_jm1)/(2*delta_t);
deedt(i,j)=(ee_i_jp1-ee_i_jm1)/(2*delta_t);
% Substituting into equations as they appear on the paper
% Equation 10
dddt(i,j) = B.*((dd_i_jp1-2.*dd(i,j)+dd_i_jm1)/(delta_x.^2)); % Equation 10
dedt(i,j) = C.*d2edx2(i,j); % Equation 11
daadt(i,j) + dffdt(i,j) + dgdt(i,j)== A.*d2aadx2(i,j) + E.*d2ffdx2(i,j)+ TT.*d2gdx2(i,j); % Equation 12
daadt(i,j) = A.*d2aadx2(i,j) - (H.*I.*L - H/M.*J.*K)*h(i,j).*P; % Equation 13
UU.*B.*dddx(i,j)+ VV.*C.*dedx(i,j) + WW.*AAA.*diidx(i,j)...
- XX.*DD.*djjdx(i,j) - YY.*E.*dffdx(i,j) - ZZ.* TT.*dgdx(i,j) == 0; % Equation 19
%====================================================================================================================================================================
elseif j>50 & j<= 98
% previous time point (i-1) and forward spatial point (j+1)
aa_im1_jp1=aa(i-1,j+1);
b_im1_jp1=b(i-1,j+1);
c_im1_jp1=c(i-1,j+1);
d_im1_jp1=d(i-1,j+1);
e_im1_jp1=e(i-1,j+1);
ff_im1_jp1=ff(i-1,j+1);
g_im1_jp1=g(i-1,j+1);
h_im1_jp1=h(i-1,j+1);
ii_im1_jp1=ii(i-1,j+1);
jj_im1_jp1=jj(i-1,j+1);
k_im1_jp1=k(i-1,j+1);
ll_im1_jp1=ll(i-1,j+1);
m_im1_jp1=m(i-1,j+1);
nn_im1_jp1=nn(i-1,j+1);
o_im1_jp1=o(i-1,j+1);
p_im1_jp1=p(i-1,j+1);
q_im1_jp1=q(i-1,j+1);
rr_im1_jp1=rr(i-1,j+1);
s_im1_jp1=s(i-1,j+1);
tt_im1_jp1=tt(i-1,j+1);
u_im1_jp1=u(i-1,j+1);
v_im1_jp1=v(i-1,j+1);
w_im1_jp1=w(i-1,j+1);
xx_im1_jp1=xx(i-1,j+1);
yy_im1_jp1=yy(i-1,j+1);
zz_im1_jp1=zz(i-1,j+1);
aaa_im1_jp1=aaa(i-1,j+1);
bb_im1_jp1=bb(i-1,j+1);
cc_im1_jp1=cc(i-1,j+1);
dd_im1_jp1=dd(i-1,j+1);
ee_im1_jp1=ee(i-1,j+1);
% previous time point (i-1) and current spatial point (j)
aa_im1_j=aa(i-1,j);
b_im1_j=b(i-1,j);
c_im1_j=c(i-1,j);
d_im1_j=d(i-1,j);
e_im1_j=e(i-1,j);
ff_im1_j=ff(i-1,j);
g_im1_j=g(i-1,j);
h_im1_j=h(i-1,j);
ii_im1_j=ii(i-1,j);
jj_im1_j=jj(i-1,j);
k_im1_j=k(i-1,j);
ll_im1_j=ll(i-1,j);
m_im1_j=m(i-1,j);
nn_im1_j=nn(i-1,j);
o_im1_j=o(i-1,j);
p_im1_j=p(i-1,j);
q_im1_j=q(i-1,j);
rr_im1_j=rr(i-1,j);
s_im1_j=s(i-1,j);
tt_im1_j=tt(i-1,j);
u_im1_j=u(i-1,j);
v_im1_j=v(i-1,j);
w_im1_j=w(i-1,j);
xx_im1_j=xx(i-1,j);
yy_im1_j=yy(i-1,j);
zz_im1_j=zz(i-1,j);
aaa_im1_j=aaa(i-1,j);
bb_im1_j=bb(i-1,j);
cc_im1_j=cc(i-1,j);
dd_im1_j=dd(i-1,j);
ee_im1_j=ee(i-1,j);
% previous time point (i-1) and previous spatial point (j-1)
aa_im1_jm1=aa(i-1,j-1);
b_im1_jm1=b(i-1,j-1);
c_im1_jm1=c(i-1,j-1);
d_im1_jm1=d(i-1,j-1);
e_im1_jm1=e(i-1,j-1);
ff_im1_jm1=ff(i-1,j-1);
g_im1_jm1=g(i-1,j-1);
h_im1_jm1=h(i-1,j-1);
ii_im1_jm1=ii(i-1,j-1);
jj_im1_jm1=jj(i-1,j-1);
k_im1_jm1=k(i-1,j-1);
ll_im1_jm1=ll(i-1,j-1);
m_im1_jm1=m(i-1,j-1);
nn_im1_jm1=nn(i-1,j-1);
o_im1_jm1=o(i-1,j-1);
p_im1_jm1=p(i-1,j-1);
q_im1_jm1=q(i-1,j-1);
rr_im1_jm1=rr(i-1,j-1);
s_im1_jm1=s(i-1,j-1);
tt_im1_jm1=tt(i-1,j-1);
u_im1_jm1=u(i-1,j-1);
v_im1_jm1=v(i-1,j-1);
w_im1_jm1=w(i-1,j-1);
xx_im1_jm1=xx(i-1,j-1);
yy_im1_jm1=yy(i-1,j-1);
zz_im1_jm1=zz(i-1,j-1);
aaa_im1_jm1=aaa(i-1,j-1);
bb_im1_jm1=bb(i-1,j-1);
cc_im1_jm1=cc(i-1,j-1);
dd_im1_jm1=dd(i-1,j-1);
ee_im1_jm1=ee(i-1,j-1);
% Substituting equations as they appear in the paper
for n = vinitial:vfinal
n(v,i) = (1/3) + g1*exp(-(1/4)*LambertW((-a6*exp((4*sqrt(a2)*alpha)/(2*a2)))) + ...
((6*g1^2)*(exp(-(1/4)*LambertW((-a6*exp(4*sqrt(a2)*alpha)/(2*a2))+sqrt(a2)*alpha)^2))/(18*v + 1))); % Equation 26
dkdt(i,j) = -a.*B.*dkdx(i,j) + integral((II.*kMT.*Ap.*(ll(i,j) - k(i,j))).*n(v),vfinal,vinitial); % Equation 20
dmdt(i,j) = -a.*C.*dmdx(i,j) + integral((II.*kMT.*Ap.*(nn(i,j) - m(i,j))).*n(v),vfinal,vinitial); % Equation 28
(dodt(i,j) + Wdt(i,j) + dqdt(i,j)) ==( -a.*A.*dodx(i,j) - a.*E.*Wdx(i,j) - a.*TT.*dqdx(i,j) - ...
((integral((II.*kMT.*Ap.*(ll(i,j) - k(i,j))).*n(v),vfinal,vinitial))+(integral((II.*kMT.*Ap.*(nn(i,j) - ...
m(i,j))).*n(v),vfinal,vinitial)))); % Equation 30
dodt(i,j) = a.*A.*d2odx2(i,j)-(H.*I.*L - H/M.*J.*K)*rr(i,j).*P; %Equation 31
dkdt(i,j)+ dmdt(i,j)+ dsdt(i,j) - dttdt(i,j) - Wdt(i,j) + dqdt(i,j) == a.*B.*dkdx(i,j) + ...
a.*C.*dmdx(i,j) + a.*AAA.*dsdx(i,j)- a.*DD.*dttdx(i,j) - a.*E.*Wdx(i,j)- a.*TT.*dqdx(i,j); %Equation 31a
end
%==================================================================================================================================================================
elseif j>100 & j<=150
% previous time point (i-1) and forward spatial point (j+1)
aa_im1_jp1=aa(i-1,j+1);
b_im1_jp1=b(i-1,j+1);
c_im1_jp1=c(i-1,j+1);
d_im1_jp1=d(i-1,j+1);
e_im1_jp1=e(i-1,j+1);
ff_im1_jp1=ff(i-1,j+1);
g_im1_jp1=g(i-1,j+1);
h_im1_jp1=h(i-1,j+1);
ii_im1_jp1=ii(i-1,j+1);
jj_im1_jp1=jj(i-1,j+1);
k_im1_jp1=k(i-1,j+1);
ll_im1_jp1=ll(i-1,j+1);
m_im1_jp1=m(i-1,j+1);
nn_im1_jp1=nn(i-1,j+1);
o_im1_jp1=o(i-1,j+1);
p_im1_jp1=p(i-1,j+1);
q_im1_jp1=q(i-1,j+1);
rr_im1_jp1=rr(i-1,j+1);
s_im1_jp1=s(i-1,j+1);
tt_im1_jp1=tt(i-1,j+1);
u_im1_jp1=u(i-1,j+1);
v_im1_jp1=v(i-1,j+1);
w_im1_jp1=w(i-1,j+1);
xx_im1_jp1=xx(i-1,j+1);
yy_im1_jp1=yy(i-1,j+1);
zz_im1_jp1=zz(i-1,j+1);
aaa_im1_jp1=aaa(i-1,j+1);
bb_im1_jp1=bb(i-1,j+1);
cc_im1_jp1=cc(i-1,j+1);
dd_im1_jp1=dd(i-1,j+1);
ee_im1_jp1=ee(i-1,j+1);
% previous time point (i-1) and current spatial point (j)
aa_im1_j=aa(i-1,j);
b_im1_j=b(i-1,j);
c_im1_j=c(i-1,j);
d_im1_j=d(i-1,j);
e_im1_j=e(i-1,j);
ff_im1_j=ff(i-1,j);
g_im1_j=g(i-1,j);
h_im1_j=h(i-1,j);
ii_im1_j=ii(i-1,j);
jj_im1_j=jj(i-1,j);
k_im1_j=k(i-1,j);
ll_im1_j=ll(i-1,j);
m_im1_j=m(i-1,j);
nn_im1_j=nn(i-1,j);
o_im1_j=o(i-1,j);
p_im1_j=p(i-1,j);
q_im1_j=q(i-1,j);
rr_im1_j=rr(i-1,j);
s_im1_j=s(i-1,j);
tt_im1_j=tt(i-1,j);
u_im1_j=u(i-1,j);
v_im1_j=v(i-1,j);
w_im1_j=w(i-1,j);
xx_im1_j=xx(i-1,j);
yy_im1_j=yy(i-1,j);
zz_im1_j=zz(i-1,j);
aaa_im1_j=aaa(i-1,j);
bb_im1_j=bb(i-1,j);
cc_im1_j=cc(i-1,j);
dd_im1_j=dd(i-1,j);
ee_im1_j=ee(i-1,j);
% previous time point (i-1) and previous spatial point (j-1)
aa_im1_jm1=aa(i-1,j-1);
b_im1_jm1=b(i-1,j-1);
c_im1_jm1=c(i-1,j-1);
d_im1_jm1=d(i-1,j-1);
e_im1_jm1=e(i-1,j-1);
ff_im1_jm1=ff(i-1,j-1);
g_im1_jm1=g(i-1,j-1);
h_im1_jm1=h(i-1,j-1);
ii_im1_jm1=ii(i-1,j-1);
jj_im1_jm1=jj(i-1,j-1);
k_im1_jm1=k(i-1,j-1);
ll_im1_jm1=ll(i-1,j-1);
m_im1_jm1=m(i-1,j-1);
nn_im1_jm1=nn(i-1,j-1);
o_im1_jm1=o(i-1,j-1);
p_im1_jm1=p(i-1,j-1);
q_im1_jm1=q(i-1,j-1);
rr_im1_jm1=rr(i-1,j-1);
s_im1_jm1=s(i-1,j-1);
tt_im1_jm1=tt(i-1,j-1);
u_im1_jm1=u(i-1,j-1);
v_im1_jm1=v(i-1,j-1);
w_im1_jm1=w(i-1,j-1);
xx_im1_jm1=xx(i-1,j-1);
yy_im1_jm1=yy(i-1,j-1);
zz_im1_jm1=zz(i-1,j-1);
aaa_im1_jm1=aaa(i-1,j-1);
bb_im1_jm1=bb(i-1,j-1);
cc_im1_jm1=cc(i-1,j-1);
dd_im1_jm1=dd(i-1,j-1);
ee_im1_jm1=ee(i-1,j-1);
% Substituting equations as they appear in the paper
4/W .* A.*d2wdx2(i,j) - rr(i,j).*P(H.*I.*x.*L.*x - kb.*J.*x.*K.*x)==0; %Equation 36
daaadx(i,150) = 0; %Equation 37
w(i,100) = o(i,j) %Equation 38
((u(i,j)/(u(i,j)+v(i,j)+w(i,j)+xx(i,j)+yy(i,j)+dd(i,j)+ee(i,j))) + ...
(v(i,j)/(u(i,j)+v(i,j)+w(i,j)+xx(i,j)+yy(i,j)+dd(i,j)+ee(i,j)))+ ...
(w(i,j)/(u(i,j)+v(i,j)+w(i,j)+xx(i,j)+yy(i,j)+dd(i,j)+ee(i,j)))+ ...
(xx(i,j)/(u(i,j)+v(i,j)+w(i,j)+xx(i,j)+yy(i,j)+dd(i,j)+ee(i,j)))+ ...
(yy(i,j)/(u(i,j)+v(i,j)+w(i,j)+xx(i,j)+yy(i,j)+dd(i,j)+ee(i,j)))+ ...
(dd(i,j)/(u(i,j)+v(i,j)+w(i,j)+xx(i,j)+yy(i,j)+dd(i,j)+ee(i,j)))+ ...
(ee(i,j)/(u(i,j)+v(i,j)+w(i,j)+xx(i,j)+yy(i,j)+dd(i,j)+ee(i,j)))) == 1; %Equation 40
%================================================================================================================================================
else j==97
% previous time point (i-1) and forward spatial point (j+1)
aa_im1_jp1=aa(i-1,j+1);
b_im1_jp1=b(i-1,j+1);
c_im1_jp1=c(i-1,j+1);
d_im1_jp1=d(i-1,j+1);
e_im1_jp1=e(i-1,j+1);
ff_im1_jp1=ff(i-1,j+1);
g_im1_jp1=g(i-1,j+1);
h_im1_jp1=h(i-1,j+1);
ii_im1_jp1=ii(i-1,j+1);
jj_im1_jp1=jj(i-1,j+1);
k_im1_jp1=k(i-1,j+1);
ll_im1_jp1=ll(i-1,j+1);
m_im1_jp1=m(i-1,j+1);
nn_im1_jp1=nn(i-1,j+1);
o_im1_jp1=o(i-1,j+1);
p_im1_jp1=p(i-1,j+1);
q_im1_jp1=q(i-1,j+1);
rr_im1_jp1=rr(i-1,j+1);
s_im1_jp1=s(i-1,j+1);
tt_im1_jp1=tt(i-1,j+1);
u_im1_jp1=u(i-1,j+1);
v_im1_jp1=v(i-1,j+1);
w_im1_jp1=w(i-1,j+1);
xx_im1_jp1=xx(i-1,j+1);
yy_im1_jp1=yy(i-1,j+1);
zz_im1_jp1=zz(i-1,j+1);
aaa_im1_jp1=aaa(i-1,j+1);
bb_im1_jp1=bb(i-1,j+1);
cc_im1_jp1=cc(i-1,j+1);
dd_im1_jp1=dd(i-1,j+1);
ee_im1_jp1=ee(i-1,j+1);
% previous time point (i-1) and current spatial point (j)
aa_im1_j=aa(i-1,j);
b_im1_j=b(i-1,j);
c_im1_j=c(i-1,j);
d_im1_j=d(i-1,j);
e_im1_j=e(i-1,j);
ff_im1_j=ff(i-1,j);
g_im1_j=g(i-1,j);
h_im1_j=h(i-1,j);
ii_im1_j=ii(i-1,j);
jj_im1_j=jj(i-1,j);
k_im1_j=k(i-1,j);
ll_im1_j=ll(i-1,j);
m_im1_j=m(i-1,j);
nn_im1_j=nn(i-1,j);
o_im1_j=o(i-1,j);
p_im1_j=p(i-1,j);
q_im1_j=q(i-1,j);
rr_im1_j=rr(i-1,j);
s_im1_j=s(i-1,j);
tt_im1_j=tt(i-1,j);
u_im1_j=u(i-1,j);
v_im1_j=v(i-1,j);
w_im1_j=w(i-1,j);
xx_im1_j=xx(i-1,j);
yy_im1_j=yy(i-1,j);
zz_im1_j=zz(i-1,j);
aaa_im1_j=aaa(i-1,j);
bb_im1_j=bb(i-1,j);
cc_im1_j=cc(i-1,j);
dd_im1_j=dd(i-1,j);
ee_im1_j=ee(i-1,j);
% previous time point (i-1) and previous spatial point (j-1)
aa_im1_jm1=aa(i-1,j-1);
b_im1_jm1=b(i-1,j-1);
c_im1_jm1=c(i-1,j-1);
d_im1_jm1=d(i-1,j-1);
e_im1_jm1=e(i-1,j-1);
ff_im1_jm1=ff(i-1,j-1);
g_im1_jm1=g(i-1,j-1);
h_im1_jm1=h(i-1,j-1);
ii_im1_jm1=ii(i-1,j-1);
jj_im1_jm1=jj(i-1,j-1);
k_im1_jm1=k(i-1,j-1);
ll_im1_jm1=ll(i-1,j-1);
m_im1_jm1=m(i-1,j-1);
nn_im1_jm1=nn(i-1,j-1);
o_im1_jm1=o(i-1,j-1);
p_im1_jm1=p(i-1,j-1);
q_im1_jm1=q(i-1,j-1);
rr_im1_jm1=rr(i-1,j-1);
s_im1_jm1=s(i-1,j-1);
tt_im1_jm1=tt(i-1,j-1);
u_im1_jm1=u(i-1,j-1);
v_im1_jm1=v(i-1,j-1);
w_im1_jm1=w(i-1,j-1);
xx_im1_jm1=xx(i-1,j-1);
yy_im1_jm1=yy(i-1,j-1);
zz_im1_jm1=zz(i-1,j-1);
aaa_im1_jm1=aaa(i-1,j-1);
bb_im1_jm1=bb(i-1,j-1);
cc_im1_jm1=cc(i-1,j-1);
dd_im1_jm1=dd(i-1,j-1);
ee_im1_jm1=ee(i-1,j-1);
((A.*(w(i,j).*x - o(i,j)) + E.*(xx(i,j).*x - p(i,j)) + TT.*(yy(i,j).*x - q(i,j)))) == ...
((B.*(ll(i,j) - k(i,j))+ C.*(nn(i,j) - m(i,j))).*(1-x)); % Equation 42
end
end
end
What can I do?
Regards
Dursman
  2 Comments
KSSV
KSSV on 19 Nov 2020
You have made your code very complex. It is tiresome to check. But note that this kind of error comes when you try to extract the more number of elements than present.
A = rand(1,5) ; % five elements are present
A(2) % no error
A(5) % no error
A(6) % error, becuase there is no sixth element

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on Geometry and Mesh in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!