I want to differentiate NAA with respect to U(L,J+1), how to do it.
%———————SPLINE IN COMPRESSION————————————————
% Numerical Solution of Burger’s Huxley equation
%E*Uxx=Ut+ALPHA*UUx+BETA*(U^3+DELTA*U^2+GAMMA*U)
% , U(x,t)=GAMMA*[1+tanh(a1*(x-a2*t))]/2, t>=0
%————————————————————————–
clear all; % clear all variables in memory
LAM=1.6;
%LL=10;
%JJ=1000; % Solution at t=1
LL=8;
JJ=40;
H=1/LL;
K=LAM*H*H;
K=0.01;
%LLL=(LL/10)+9;
H= (2*pi)/100;
SGM=1.0;
EPS=1.0;
ALPHA=1.0;
BETA=0.0;
GAMMA=0.0;
DELTA=1.0;
a= (-1/4)*(1-SGM+(SGM)^2);
b= -SGM/(3*(1+SGM));
c= (-1/(2*SGM*(1+SGM)))*(1-SGM+(SGM)^2);
%a1=GAMMA*(-ALPHA+sqrt((ALPHA^2)+8*BETA))/8;
%a2=(ALPHA*GAMMA/2)+((2-GAMMA)*(ALPHA+sqrt((ALPHA^2)+8*BETA))/4);
%a1=(-ALPHA-sqrt((ALPHA^2)+8*BETA))/8;
%a2=(ALPHA/2)-((1-2*GAMMA)*(ALPHA-sqrt((ALPHA^2)+8*BETA))/4);
% LLL=(0.5*LL)+1;
%tt=(8.986818916*pi)/180;
%tt=8.986818916;
% %Spline in Tension
% tt=0.001;
% AA=(SGM*((tt/sin(tt))-cos(tt)))/(2*tt^2);
% BB1=(SGM*(-(tt*cos(tt)/sin(tt))+1))/(2*tt^2);
% BB2=(-(tt*cos(tt)/sin(tt))+1)/(2*tt^2);
% CC=((tt/sin(tt))-cos(tt))/(2*tt^2);
% %spline in compression
tt=(8.986818916*pi)/180;
AA=(SGM*((tt/sin(tt))-cos(tt)))/(2*tt^2);
BB1=(SGM*(-(tt*cos(tt)/sin(tt))+1))/(2*tt^2);
BB2=(-(tt*cos(tt)/sin(tt))+1)/(2*tt^2);
CC=((tt/sin(tt))-cos(tt))/(2*tt^2);
%Cubic Spline
%AA=SGM/3;
%BB1=SGM/6;
%BB2=1/6;
%CC=1/3;
disp(‘Numerical Solution of 1D Parabolic Equation’);
disp(‘ ‘);
fprintf(‘LL=%4.2f, K=%6.4f, LAM=%6.2f n’,LL,K,LAM);
disp(‘ ‘);
% Dimensions
X=zeros(1,LL+1);
T=zeros(1,JJ+1);
TT=zeros(1,JJ+1);
U=zeros(LL+1,JJ+1);
EU=zeros(LL+1,JJ+1);
OLDU=zeros(LL+1,JJ+1);
H=zeros(1,LL);
if SGM==1
S=1;
for i=1:LL-1
S=S+SGM^i;
end;
% H(1)=2/S;
H(1)=1/S;
else
% H(1)=2*((1-SGM)/(1-(SGM^LL)));
H(1)=((1-SGM)/(1-(SGM^LL)));
end
% Formation of Variable step lengths
for L=1:LL-1
H(L+1)=SGM*H(L);
end
% Formation of Nodal Points
% X(1)=-1.0;
X(1)=0.0;
for L=2:LL
for J=1:JJ+1
X(L)=X(L-1)+H(L-1);
T(J)=(J-1)*K;
end
end
X(LL+1)=1.0;
for L=1:LL+1
for J=1:JJ+1
% EU(L,J)=GAMMA*(1+tanh(a1*(X(L)-a2*T(J))))/2;
% EU(L,J)=(1+tanh(a1*(X(L)-a2*T(J))))/2;
EU(L,J)= exp(-T(J))*sin(X(L));
end
end
% Initial conditions
for L=1:LL+1
%U(L,1)=EU(L,1);
U(L,1)=cos(pi*X(L));
end
% Boundary conditions
for J=2:JJ+1
%U(0,J)==EU(0,J);
U(1,J)=exp(-EPS*(pi)^2*T(J));
U(2,J)=-exp(-EPS*(pi)^2*T(J));
%U(1,J)=EU(1,J)==-exp(-EPS*(pi)^2*t);;
U(LL+1,J)=EU(LL+1,J);
end
% Solution at advanced time level
tic % CPU time starts
for J=1:JJ
% Initial & First Approx
for L=2:LL
OLDU(L,J+1)=EU(L,J+1);
U(L,J+1)=OLDU(L,J+1);
end
for KK=1:1200 %Iteration starts
% Solution at Internal grid points at each time level
for L=2:LL
M1 = U(L+1,J+1);
M2 = U(L+1,J);
M3 = U(L,J+1);
M4 = U(L,J);
M5 = U(L-1,J+1);
M6 = U(L-1,J);
N1 = (1/(2*K))*(M5+M3-M6-M4);
N2 = (M5+M3+M6+M4)*(M5+M3+M6+M4-2)*(M5+M3+M6+M4-2*GAMMA);
N3 = (M5+M3++M6+M4);
N4 = (1/(2*SGM*H(L)))*(M5+M6-M3-M4);
N5 = (1/K)*(M3-M4)+(BETA/8)*(M3+M4)*(M3+M4-2)*(M3+M4-2*GAMMA)+(ALPHA/2)*(M3+M4)*(1/(2*SGM*H(L)*(SGM+1)))*(M5+M6-(1-(SGM)^2)*(M3+M4)-((SGM)^2)*(M1+M2));
N6 = (1/(2*K))*(M5+M3-M6-M4)+(ALPHA/8)*(M5+M6+M3+M4)*(M5+M6+M3+M4)+(BETA/64)*(M5+M6+M3+M4)*(M5+M6+M3+M4-4)*(M5+M6+M3+M4-4*GAMMA);
N7 = (1/(2*K))*(M1+M3-M2-M4);
N8 = (M1+M3+M2+M4)*(M1+M3+M2+M4-2)*(M1+M3+M2+M4-2*GAMMA);
N9 = (M1+M3++M2+M4);
N10 = (1/(2*SGM*H(L)))*(M1+M2-M3-M4);
N11 = (1/K)*(M3-M4)+(BETA/8)*(M3+M4)*(M3+M4-2)*(M3+M4-2*GAMMA)+(ALPHA/2)*(M3+M4)*(1/(2*SGM*H(L)*(SGM+1)))*(M1+M2-(1-(SGM)^2)*(M3+M4)-((SGM)^2)*(M5+M6));
N12 = (1/(2*K))*(M1+M3-M2-M4)+(ALPHA/8)*(M1+M2+M3+M4)*(M1+M2+M3+M4)+(BETA/64)*(M1+M2+M3+M4)*(M1+M2+M3+M4-4)*(M1+M2+M3+M4-4*GAMMA);
N13 = (1/K)*(M3-M4);
N14 = (1/K)*(M1-M2)-(1+SGM)*N13+(SGM/K)*(M5-M6);
N15= (M3+M4)/(2);
N16 = M1+M2+(1+SGM)*(M3+M4)+SGM*(M5+M6);
N17 = M1+M2-(1-(SGM^2))*(M3+M4)-(SGM)^2*(M5+M6);
N18 = (1/(2*K))*(M1+M3-M2-M4);
N19 =N9^2;
N21 = M5+M3-M6-M4;
N22 = N3^2;
O1 = N1+(BETA/16)*(N2)+(ALPHA/4)*(N3)*(N4+(H(L)/4)*(2*BB2*N5-CC*N6));
O2 = N7+(BETA/16)*(N8)+(ALPHA/4)*(N9)*(N10+(H(L)/4)*(2*BB1*N11-AA*N12));
O3 = N13+c*N14+ALPHA*(N15+(a/SGM*(1+SGM))*N16)*(1/2*H(L)*SGM*(1+SGM))*N17 + b*H(L)*[N18+(ALPHA/8)*N19 + (BETA/16)*N9*(N9-4)*(N9-4*GAMMA)-(N21+(ALPHA/8)*(N3)^2)+(BETA/16)*(N3)*(N3-4)*(N3-4*GAMMA)]+BETA*[(N15+(a/(SGM*(SGM+1)))*N16)*(N15+(a/(2*SGM*(SGM+1)))*N16-1)*(N15+(a/(2*SGM*(SGM+1)))*N16-GAMMA)];
NAA = (EPS/2)*(M1+M2-2*M3-2*M4+M5+M6) – ((H.^2)/3)* (O1+O2+O3);%———————SPLINE IN COMPRESSION————————————————
% Numerical Solution of Burger’s Huxley equation
%E*Uxx=Ut+ALPHA*UUx+BETA*(U^3+DELTA*U^2+GAMMA*U)
% , U(x,t)=GAMMA*[1+tanh(a1*(x-a2*t))]/2, t>=0
%————————————————————————–
clear all; % clear all variables in memory
LAM=1.6;
%LL=10;
%JJ=1000; % Solution at t=1
LL=8;
JJ=40;
H=1/LL;
K=LAM*H*H;
K=0.01;
%LLL=(LL/10)+9;
H= (2*pi)/100;
SGM=1.0;
EPS=1.0;
ALPHA=1.0;
BETA=0.0;
GAMMA=0.0;
DELTA=1.0;
a= (-1/4)*(1-SGM+(SGM)^2);
b= -SGM/(3*(1+SGM));
c= (-1/(2*SGM*(1+SGM)))*(1-SGM+(SGM)^2);
%a1=GAMMA*(-ALPHA+sqrt((ALPHA^2)+8*BETA))/8;
%a2=(ALPHA*GAMMA/2)+((2-GAMMA)*(ALPHA+sqrt((ALPHA^2)+8*BETA))/4);
%a1=(-ALPHA-sqrt((ALPHA^2)+8*BETA))/8;
%a2=(ALPHA/2)-((1-2*GAMMA)*(ALPHA-sqrt((ALPHA^2)+8*BETA))/4);
% LLL=(0.5*LL)+1;
%tt=(8.986818916*pi)/180;
%tt=8.986818916;
% %Spline in Tension
% tt=0.001;
% AA=(SGM*((tt/sin(tt))-cos(tt)))/(2*tt^2);
% BB1=(SGM*(-(tt*cos(tt)/sin(tt))+1))/(2*tt^2);
% BB2=(-(tt*cos(tt)/sin(tt))+1)/(2*tt^2);
% CC=((tt/sin(tt))-cos(tt))/(2*tt^2);
% %spline in compression
tt=(8.986818916*pi)/180;
AA=(SGM*((tt/sin(tt))-cos(tt)))/(2*tt^2);
BB1=(SGM*(-(tt*cos(tt)/sin(tt))+1))/(2*tt^2);
BB2=(-(tt*cos(tt)/sin(tt))+1)/(2*tt^2);
CC=((tt/sin(tt))-cos(tt))/(2*tt^2);
%Cubic Spline
%AA=SGM/3;
%BB1=SGM/6;
%BB2=1/6;
%CC=1/3;
disp(‘Numerical Solution of 1D Parabolic Equation’);
disp(‘ ‘);
fprintf(‘LL=%4.2f, K=%6.4f, LAM=%6.2f n’,LL,K,LAM);
disp(‘ ‘);
% Dimensions
X=zeros(1,LL+1);
T=zeros(1,JJ+1);
TT=zeros(1,JJ+1);
U=zeros(LL+1,JJ+1);
EU=zeros(LL+1,JJ+1);
OLDU=zeros(LL+1,JJ+1);
H=zeros(1,LL);
if SGM==1
S=1;
for i=1:LL-1
S=S+SGM^i;
end;
% H(1)=2/S;
H(1)=1/S;
else
% H(1)=2*((1-SGM)/(1-(SGM^LL)));
H(1)=((1-SGM)/(1-(SGM^LL)));
end
% Formation of Variable step lengths
for L=1:LL-1
H(L+1)=SGM*H(L);
end
% Formation of Nodal Points
% X(1)=-1.0;
X(1)=0.0;
for L=2:LL
for J=1:JJ+1
X(L)=X(L-1)+H(L-1);
T(J)=(J-1)*K;
end
end
X(LL+1)=1.0;
for L=1:LL+1
for J=1:JJ+1
% EU(L,J)=GAMMA*(1+tanh(a1*(X(L)-a2*T(J))))/2;
% EU(L,J)=(1+tanh(a1*(X(L)-a2*T(J))))/2;
EU(L,J)= exp(-T(J))*sin(X(L));
end
end
% Initial conditions
for L=1:LL+1
%U(L,1)=EU(L,1);
U(L,1)=cos(pi*X(L));
end
% Boundary conditions
for J=2:JJ+1
%U(0,J)==EU(0,J);
U(1,J)=exp(-EPS*(pi)^2*T(J));
U(2,J)=-exp(-EPS*(pi)^2*T(J));
%U(1,J)=EU(1,J)==-exp(-EPS*(pi)^2*t);;
U(LL+1,J)=EU(LL+1,J);
end
% Solution at advanced time level
tic % CPU time starts
for J=1:JJ
% Initial & First Approx
for L=2:LL
OLDU(L,J+1)=EU(L,J+1);
U(L,J+1)=OLDU(L,J+1);
end
for KK=1:1200 %Iteration starts
% Solution at Internal grid points at each time level
for L=2:LL
M1 = U(L+1,J+1);
M2 = U(L+1,J);
M3 = U(L,J+1);
M4 = U(L,J);
M5 = U(L-1,J+1);
M6 = U(L-1,J);
N1 = (1/(2*K))*(M5+M3-M6-M4);
N2 = (M5+M3+M6+M4)*(M5+M3+M6+M4-2)*(M5+M3+M6+M4-2*GAMMA);
N3 = (M5+M3++M6+M4);
N4 = (1/(2*SGM*H(L)))*(M5+M6-M3-M4);
N5 = (1/K)*(M3-M4)+(BETA/8)*(M3+M4)*(M3+M4-2)*(M3+M4-2*GAMMA)+(ALPHA/2)*(M3+M4)*(1/(2*SGM*H(L)*(SGM+1)))*(M5+M6-(1-(SGM)^2)*(M3+M4)-((SGM)^2)*(M1+M2));
N6 = (1/(2*K))*(M5+M3-M6-M4)+(ALPHA/8)*(M5+M6+M3+M4)*(M5+M6+M3+M4)+(BETA/64)*(M5+M6+M3+M4)*(M5+M6+M3+M4-4)*(M5+M6+M3+M4-4*GAMMA);
N7 = (1/(2*K))*(M1+M3-M2-M4);
N8 = (M1+M3+M2+M4)*(M1+M3+M2+M4-2)*(M1+M3+M2+M4-2*GAMMA);
N9 = (M1+M3++M2+M4);
N10 = (1/(2*SGM*H(L)))*(M1+M2-M3-M4);
N11 = (1/K)*(M3-M4)+(BETA/8)*(M3+M4)*(M3+M4-2)*(M3+M4-2*GAMMA)+(ALPHA/2)*(M3+M4)*(1/(2*SGM*H(L)*(SGM+1)))*(M1+M2-(1-(SGM)^2)*(M3+M4)-((SGM)^2)*(M5+M6));
N12 = (1/(2*K))*(M1+M3-M2-M4)+(ALPHA/8)*(M1+M2+M3+M4)*(M1+M2+M3+M4)+(BETA/64)*(M1+M2+M3+M4)*(M1+M2+M3+M4-4)*(M1+M2+M3+M4-4*GAMMA);
N13 = (1/K)*(M3-M4);
N14 = (1/K)*(M1-M2)-(1+SGM)*N13+(SGM/K)*(M5-M6);
N15= (M3+M4)/(2);
N16 = M1+M2+(1+SGM)*(M3+M4)+SGM*(M5+M6);
N17 = M1+M2-(1-(SGM^2))*(M3+M4)-(SGM)^2*(M5+M6);
N18 = (1/(2*K))*(M1+M3-M2-M4);
N19 =N9^2;
N21 = M5+M3-M6-M4;
N22 = N3^2;
O1 = N1+(BETA/16)*(N2)+(ALPHA/4)*(N3)*(N4+(H(L)/4)*(2*BB2*N5-CC*N6));
O2 = N7+(BETA/16)*(N8)+(ALPHA/4)*(N9)*(N10+(H(L)/4)*(2*BB1*N11-AA*N12));
O3 = N13+c*N14+ALPHA*(N15+(a/SGM*(1+SGM))*N16)*(1/2*H(L)*SGM*(1+SGM))*N17 + b*H(L)*[N18+(ALPHA/8)*N19 + (BETA/16)*N9*(N9-4)*(N9-4*GAMMA)-(N21+(ALPHA/8)*(N3)^2)+(BETA/16)*(N3)*(N3-4)*(N3-4*GAMMA)]+BETA*[(N15+(a/(SGM*(SGM+1)))*N16)*(N15+(a/(2*SGM*(SGM+1)))*N16-1)*(N15+(a/(2*SGM*(SGM+1)))*N16-GAMMA)];
NAA = (EPS/2)*(M1+M2-2*M3-2*M4+M5+M6) – ((H.^2)/3)* (O1+O2+O3); %———————SPLINE IN COMPRESSION————————————————
% Numerical Solution of Burger’s Huxley equation
%E*Uxx=Ut+ALPHA*UUx+BETA*(U^3+DELTA*U^2+GAMMA*U)
% , U(x,t)=GAMMA*[1+tanh(a1*(x-a2*t))]/2, t>=0
%————————————————————————–
clear all; % clear all variables in memory
LAM=1.6;
%LL=10;
%JJ=1000; % Solution at t=1
LL=8;
JJ=40;
H=1/LL;
K=LAM*H*H;
K=0.01;
%LLL=(LL/10)+9;
H= (2*pi)/100;
SGM=1.0;
EPS=1.0;
ALPHA=1.0;
BETA=0.0;
GAMMA=0.0;
DELTA=1.0;
a= (-1/4)*(1-SGM+(SGM)^2);
b= -SGM/(3*(1+SGM));
c= (-1/(2*SGM*(1+SGM)))*(1-SGM+(SGM)^2);
%a1=GAMMA*(-ALPHA+sqrt((ALPHA^2)+8*BETA))/8;
%a2=(ALPHA*GAMMA/2)+((2-GAMMA)*(ALPHA+sqrt((ALPHA^2)+8*BETA))/4);
%a1=(-ALPHA-sqrt((ALPHA^2)+8*BETA))/8;
%a2=(ALPHA/2)-((1-2*GAMMA)*(ALPHA-sqrt((ALPHA^2)+8*BETA))/4);
% LLL=(0.5*LL)+1;
%tt=(8.986818916*pi)/180;
%tt=8.986818916;
% %Spline in Tension
% tt=0.001;
% AA=(SGM*((tt/sin(tt))-cos(tt)))/(2*tt^2);
% BB1=(SGM*(-(tt*cos(tt)/sin(tt))+1))/(2*tt^2);
% BB2=(-(tt*cos(tt)/sin(tt))+1)/(2*tt^2);
% CC=((tt/sin(tt))-cos(tt))/(2*tt^2);
% %spline in compression
tt=(8.986818916*pi)/180;
AA=(SGM*((tt/sin(tt))-cos(tt)))/(2*tt^2);
BB1=(SGM*(-(tt*cos(tt)/sin(tt))+1))/(2*tt^2);
BB2=(-(tt*cos(tt)/sin(tt))+1)/(2*tt^2);
CC=((tt/sin(tt))-cos(tt))/(2*tt^2);
%Cubic Spline
%AA=SGM/3;
%BB1=SGM/6;
%BB2=1/6;
%CC=1/3;
disp(‘Numerical Solution of 1D Parabolic Equation’);
disp(‘ ‘);
fprintf(‘LL=%4.2f, K=%6.4f, LAM=%6.2f n’,LL,K,LAM);
disp(‘ ‘);
% Dimensions
X=zeros(1,LL+1);
T=zeros(1,JJ+1);
TT=zeros(1,JJ+1);
U=zeros(LL+1,JJ+1);
EU=zeros(LL+1,JJ+1);
OLDU=zeros(LL+1,JJ+1);
H=zeros(1,LL);
if SGM==1
S=1;
for i=1:LL-1
S=S+SGM^i;
end;
% H(1)=2/S;
H(1)=1/S;
else
% H(1)=2*((1-SGM)/(1-(SGM^LL)));
H(1)=((1-SGM)/(1-(SGM^LL)));
end
% Formation of Variable step lengths
for L=1:LL-1
H(L+1)=SGM*H(L);
end
% Formation of Nodal Points
% X(1)=-1.0;
X(1)=0.0;
for L=2:LL
for J=1:JJ+1
X(L)=X(L-1)+H(L-1);
T(J)=(J-1)*K;
end
end
X(LL+1)=1.0;
for L=1:LL+1
for J=1:JJ+1
% EU(L,J)=GAMMA*(1+tanh(a1*(X(L)-a2*T(J))))/2;
% EU(L,J)=(1+tanh(a1*(X(L)-a2*T(J))))/2;
EU(L,J)= exp(-T(J))*sin(X(L));
end
end
% Initial conditions
for L=1:LL+1
%U(L,1)=EU(L,1);
U(L,1)=cos(pi*X(L));
end
% Boundary conditions
for J=2:JJ+1
%U(0,J)==EU(0,J);
U(1,J)=exp(-EPS*(pi)^2*T(J));
U(2,J)=-exp(-EPS*(pi)^2*T(J));
%U(1,J)=EU(1,J)==-exp(-EPS*(pi)^2*t);;
U(LL+1,J)=EU(LL+1,J);
end
% Solution at advanced time level
tic % CPU time starts
for J=1:JJ
% Initial & First Approx
for L=2:LL
OLDU(L,J+1)=EU(L,J+1);
U(L,J+1)=OLDU(L,J+1);
end
for KK=1:1200 %Iteration starts
% Solution at Internal grid points at each time level
for L=2:LL
M1 = U(L+1,J+1);
M2 = U(L+1,J);
M3 = U(L,J+1);
M4 = U(L,J);
M5 = U(L-1,J+1);
M6 = U(L-1,J);
N1 = (1/(2*K))*(M5+M3-M6-M4);
N2 = (M5+M3+M6+M4)*(M5+M3+M6+M4-2)*(M5+M3+M6+M4-2*GAMMA);
N3 = (M5+M3++M6+M4);
N4 = (1/(2*SGM*H(L)))*(M5+M6-M3-M4);
N5 = (1/K)*(M3-M4)+(BETA/8)*(M3+M4)*(M3+M4-2)*(M3+M4-2*GAMMA)+(ALPHA/2)*(M3+M4)*(1/(2*SGM*H(L)*(SGM+1)))*(M5+M6-(1-(SGM)^2)*(M3+M4)-((SGM)^2)*(M1+M2));
N6 = (1/(2*K))*(M5+M3-M6-M4)+(ALPHA/8)*(M5+M6+M3+M4)*(M5+M6+M3+M4)+(BETA/64)*(M5+M6+M3+M4)*(M5+M6+M3+M4-4)*(M5+M6+M3+M4-4*GAMMA);
N7 = (1/(2*K))*(M1+M3-M2-M4);
N8 = (M1+M3+M2+M4)*(M1+M3+M2+M4-2)*(M1+M3+M2+M4-2*GAMMA);
N9 = (M1+M3++M2+M4);
N10 = (1/(2*SGM*H(L)))*(M1+M2-M3-M4);
N11 = (1/K)*(M3-M4)+(BETA/8)*(M3+M4)*(M3+M4-2)*(M3+M4-2*GAMMA)+(ALPHA/2)*(M3+M4)*(1/(2*SGM*H(L)*(SGM+1)))*(M1+M2-(1-(SGM)^2)*(M3+M4)-((SGM)^2)*(M5+M6));
N12 = (1/(2*K))*(M1+M3-M2-M4)+(ALPHA/8)*(M1+M2+M3+M4)*(M1+M2+M3+M4)+(BETA/64)*(M1+M2+M3+M4)*(M1+M2+M3+M4-4)*(M1+M2+M3+M4-4*GAMMA);
N13 = (1/K)*(M3-M4);
N14 = (1/K)*(M1-M2)-(1+SGM)*N13+(SGM/K)*(M5-M6);
N15= (M3+M4)/(2);
N16 = M1+M2+(1+SGM)*(M3+M4)+SGM*(M5+M6);
N17 = M1+M2-(1-(SGM^2))*(M3+M4)-(SGM)^2*(M5+M6);
N18 = (1/(2*K))*(M1+M3-M2-M4);
N19 =N9^2;
N21 = M5+M3-M6-M4;
N22 = N3^2;
O1 = N1+(BETA/16)*(N2)+(ALPHA/4)*(N3)*(N4+(H(L)/4)*(2*BB2*N5-CC*N6));
O2 = N7+(BETA/16)*(N8)+(ALPHA/4)*(N9)*(N10+(H(L)/4)*(2*BB1*N11-AA*N12));
O3 = N13+c*N14+ALPHA*(N15+(a/SGM*(1+SGM))*N16)*(1/2*H(L)*SGM*(1+SGM))*N17 + b*H(L)*[N18+(ALPHA/8)*N19 + (BETA/16)*N9*(N9-4)*(N9-4*GAMMA)-(N21+(ALPHA/8)*(N3)^2)+(BETA/16)*(N3)*(N3-4)*(N3-4*GAMMA)]+BETA*[(N15+(a/(SGM*(SGM+1)))*N16)*(N15+(a/(2*SGM*(SGM+1)))*N16-1)*(N15+(a/(2*SGM*(SGM+1)))*N16-GAMMA)];
NAA = (EPS/2)*(M1+M2-2*M3-2*M4+M5+M6) – ((H.^2)/3)* (O1+O2+O3); #differentiate, #code, #pde MATLAB Answers — New Questions