After adding end to all of the functions so that this would even run, you get an dimension error (see below).
You are trying assign a 1x100 array x to a single element of X. You need to figure out what you really meant to do there.
pde2fshear_v4Perturbed_nonlinear
% solve 3-F bifurcaiton model
function pde2fshear_v4Perturbed_nonlinear
global x;
global X;
global t;
global chi0; % declare global variables
global D0;
global chi1;
global D1;
global alpha_chi;
global alpha_D;
global H0;
global S0;
global H;
global S;
global data;
global track;
global xstep;
global tstep;
global count;
global chi_growth;
global lambda_suppress;
global v_e;
global chi_ano;
global D_ano;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global B0;
global xmax;
global xmin;
%define values for constants
data.constant.critgradpressure = 1.2;
data.constant.critgraddensity = 1.0;
count = 1;
chi0 = 0.5;
D0 = 0.5;
chi1 = 5.0;
D1 = 5.0;
alpha_chi = 0.1;
alpha_D = 0.1;
H0 = 27;
S0 = 21;
track = 1;
chi_growth=20;
lambda_suppress=0.05;
chi_ano=10;
D_ano=10;
% For nonlinear model
gamma_nonlin = 50;
alpha_nonlin = 1;
D_0 = 10;
%position and time grids information
xstep = 100;
tstep = 1000;
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 50;
%tmax = 1;
m = 0; %define type of equation to solve
%Preallocate vectors for speed improvement
grad_u1 = zeros(tstep,xstep);
grad_u2 = zeros(tstep,xstep);
grad_u3 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
flowshear_p = zeros(tstep,xstep);
flowshear_n = zeros(tstep,xstep);
Q = zeros(tstep,xstep);
Q0 = zeros(tstep,xstep);
Q1 = zeros(tstep,xstep);
neo_p = zeros(tstep,xstep);
ano_p = zeros(tstep,xstep);
Gam = zeros(tstep,xstep);
Gam0 = zeros(tstep,xstep);
Gam1 = zeros(tstep,xstep);
neo_n = zeros(tstep,xstep);
ano_n = zeros(tstep,xstep);
safetyFact=zeros(tstep,xstep);
elec_diaVel=zeros(tstep,xstep);
grad_length=zeros(tstep,xstep);
mag_shear=zeros(tstep,xstep);
theta_deg=zeros(tstep,xstep);
elec_temp=zeros(tstep,xstep);
T_magFld=zeros(tstep,xstep);
P_magFld=zeros(tstep,xstep);
R0=zeros(tstep,xstep);
I_p=zeros(tstep,xstep);
jb=zeros(tstep,xstep);
B0=zeros(tstep,xstep);
%define first inner half of the plasma
%x = linspace(xmin,xmax/2,xstep/5);
x = linspace(xmin,xmax,xstep);
t = linspace(tmin,tmax,tstep);
%Add smaller grid size near plasma edge
%for i=(xstep/5)+1:xstep
% x(i) = x(i-1)+(xmax/2)/(xstep-(xstep/5));
%end
data.variable.x = x;
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
% Extract the first solution component as u1 = pressure
% second solution component as u2 = density
% third solution component as u3 = turbulence intensity
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
%grad_u = gradient(u,(x(2)-x(1)));
for j=1:tstep
for i = 1:xstep/5
if i == 1
grad_u1(j,i) = (u1(j,2)-u1(j,1))/(x(2)-x(1));
grad_u2(j,i) = (u2(j,2)-u2(j,1))/(x(2)-x(1));
grad_u3(j,i) = (u3(j,2)-u3(j,1))/(x(2)-x(1));
% grad_u4(j,i) = (u4(j,2)-u4(j,1))/(x(2)-x(1));
elseif i == xstep/5
grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
% grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
else
grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));
% grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
end
end
for i=(xstep/5)+1:xstep
if i == xstep/5+1
grad_u1(j,i) = (u1(j,i+1)-u1(j,i))/(x(i+1)-x(i));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i))/(x(i+1)-x(i));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i))/(x(i+1)-x(i));
% grad_u4(j,i) = (u4(j,i+1)-u4(j,i))/(x(i+1)-x(i));
elseif i == xstep
grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
% grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
else
grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));
% grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
end
end
end
for i=1:tstep
for j=1:xstep
v_e = -grad_u1(i,j)*grad_u2(i,j)/u2(i,j)^2; % -g_p*g_n/n^2
flowshear_p(i,j) = 1+ alpha_chi*v_e^2;
flowshear_n(i,j) = 1+ alpha_D*v_e^2;
if abs(grad_u1(i,j)) < abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
Q(i,j) = -grad_u1(i,j)*chi0;
Q0(i,j) = Q(i,j);
Q1(i,j) = 0;
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = 0;
Gam(i,j) = -grad_u2(i,j)*D0;
Gam0(i,j) = Gam(i,j);
Gam1(i,j) = 0;
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = 0;
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
elseif abs(grad_u1(i,j)) >= abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
Q0(i,j) = -grad_u1(i,j)*chi0;
Q1(i,j) = (chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = chi_ano*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
Gam(i,j) = -grad_u2(i,j)*D0;
Gam0(i,j) = Gam(i,j);
Gam1(i,j) = 0;
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = 0;
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
else
Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
Q0(i,j) = -grad_u1(i,j)*chi0;
Q1(i,j) = (chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = chi_ano*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
Gam(i,j) = -grad_u2(i,j)*(D0 + D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j));
Gam0(i,j) = -grad_u2(i,j)*D0;
Gam1(i,j) = -grad_u2(i,j)*D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j);
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = D_ano*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(i,j);
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
end
end
end
%curve_u = gradient(grad_u,(x(2)-x(1)));
for j=1:tstep
for i = 1:xstep/5
if i == 1
curve_u1(j,i) = (grad_u1(j,2)-grad_u1(j,1))/(x(2)-x(1));
curve_u2(j,i) = (grad_u2(j,2)-grad_u2(j,1))/(x(2)-x(1));
curve_u3(j,i) = (grad_u3(j,2)-grad_u3(j,1))/(x(2)-x(1));
% curve_u4(j,i) = (grad_u4(j,2)-grad_u4(j,1))/(x(2)-x(1));
elseif i == xstep/5
curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
% curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
else
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));
% curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
end
end
for i=(xstep/5)+1:xstep
if i == xstep/5+1
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i))/(x(i+1)-x(i));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i))/(x(i+1)-x(i));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i))/(x(i+1)-x(i));
% curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i))/(x(i+1)-x(i));
elseif i == xstep
curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
% curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
else
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));
% curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
end
end
end
%to save parameters and variables
data.constant.chi0 = chi0;
data.constant.D0 = D0;
data.constant.chi1 = chi1;
data.constant.D1 = D1;
data.constant.alphachi = alpha_chi;
data.constant.alphaD = alpha_D;
data.constant.H0 = H0;
data.constant.S0 = S0;
data.variable.pressure = u1;
data.variable.density = u2;
data.variable.turbulence= u3;
data.variable.intensity_diff=intensity_diff;
data.variable.drift_velFluct=drift_velFluct;
data.variable.drift_vel=drift_vel;
data.variable.gradpressure = grad_u1;
data.variable.graddensity = grad_u2;
data.variable.gradintensity = grad_u3;
data.variable.curvepressure = curve_u1;
data.variable.curvedensity = curve_u2;
data.variable.curveturbulence = curve_u3;
data.variable.x = x;
data.variable.t = t;
data.variable.Q = Q;
data.variable.Gamma = Gam;
data.variable.Q0 = Q0;
data.variable.Gamma0 = Gam0;
data.variable.neo_P = neo_p;
data.variable.neo_n = neo_n;
data.variable.Q1 = Q1;
data.variable.Gam1 = Gam1;
data.variable.ano_p = ano_p;
data.variable.ano_n = ano_n;
data.variable.heatsource = H;
data.variable.particlesource = S;
data.variable.wexb_p = flowshear_p;
data.variable.wexb_n = flowshear_n;
data.control.xgrid = xstep;
data.control.tgrid = tstep;
data.control.xmin = xmin;
data.control.xmax = xmax;
data.control.tmin = tmin;
data.control.tmax = tmax;
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u1(end,:),'.');
axis ([0 1 0 20]);
ylabel('p')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u1(end,:),'.');
axis([0 1 -60 0]);
ylabel('p\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u1(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('p\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u1(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|p\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u2(end,:),'.');
axis ([0 1 0 20]);
ylabel('n')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u2(end,:),'.');
axis([0 1 -60 0]);
ylabel('n\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u2(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('n\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u2(end,:)),Gam(end,:),'.');
%plot(nn_pre,Gam_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('\Gamma')
xlabel('|n\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u3(end,:),'.');
axis ([0 1 0 20]);
ylabel('I')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u3(end,:),'.');
axis([0 1 -60 0]);
ylabel('I\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u3(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('I\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u3(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|I\prime|')
%title('Bifurcation diagram')
grid on
end
% --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global x;
global t;
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global xstep;
global tstep;
global alpha_chi;
global alpha_D;
global chi_growth; % total growth rate
global length;
global theta_heaviside1;
global lambda_suppress;
global v_e;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global B0;
global X;
global xmax;
global xmin;
%lf FFf F Fb vfDdength = 0.01 or 1
x_count=1;
length=1;
jb0=1;
jd0=1;
x0=1;
grad_length=DuDx(2);
R0=linspace(1,5,100);
B0=1;
c = [1;1;1];
v_e = -DuDx(1)*DuDx(2)/u(2)^2; % -g_p*g_n/n^2
flowshear_p = 1+ alpha_chi*v_e^2;
flowshear_n = 1+ alpha_D*v_e^2;
X = linspace(xmin,xmax,xstep-1);
term1 = abs(DuDx(1))-data.constant.critgradpressure;
drift_velFluct = D_0*DuDx(3)^alpha_nonlin;
intensity_diff = D_0*u(3)^alpha_nonlin;
% Implementing Heaviside function for H-mode in p, n and I equations
if term1 > 0
theta_heaviside1=1;
else
theta_heaviside1=0;
end
% Determining group velocity of turbulence intensity
if x_count > xstep-1
x_count=1;
end
if x_count <= xstep-1
whos X x x_count
X(x_count)=x;
num=15;
if x_count==1
jb=jb0*DuDx(1);
jd=jd0.*((1-X(x_count)-x0).^2).^num;
j_tot=jb(x_count)+jd;
for n=1:5
I_p=trapz(x(1:n),jb(x(1:n))+jd(x(1:n)));
end
end
T_magFld=B0./R0;
P_magFld=B0.*I_p./x;
safetyFact=(x./R0).*T_magFld./P_magFld;
if x_count == 1
grad_safetyFact(x_count) = (safetyFact(2)-safetyFact(1))/(x(2)-x(1));
elseif x_count == xstep/5
grad_safetyFact(x_count) = (safetyFact(x_count)-safetyFact(x_count-1))/(x(x_count)-x(x_count-1));
else
grad_safetyFact(x_count) = (safetyFact(x_count+1)-safetyFact(x_count-1))/(x(xstep+1)-x(x_count-1));
end
mag_shear=(x/safetyFact)*grad_safetyFact;
end
elec_diaVel = -elec_temp/1.6*10^(-19)*T_magFld*grad_length;
group_vel = - (elec_diaVel*(2*grad_length/mag_shear*R0)*sin(theta_deg));
%Turbulence intensity Equation for nonlinear turbulence intensity
s3 = (chi_growth*(term1*theta_heaviside1-lambda_suppress*v_e^2)-gamma_nonlin*u(3)^alpha_nonlin)*u(3)+group_vel*u(3)-(D_0*(1+alpha_nonlin)*alpha_nonlin*u(3)^(alpha_nonlin-1)*u(3)^2);
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2;s3];
f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ;-(intensity_diff-D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin)].*DuDx; % flux term for nonlinear model
end
% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01;0.1*exp(-100*(x-1)^2)];
u0 = [0.1*(1-x^2); 0.1*(1-x^2); 0.5*exp(-100*(x-1)^2)];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = [0; 0; 0];
ql = [1; 1; 1];
pr = [ur(1); ur(2);0];
qr = [0.01; 0.1; 1];
end
%---------------------------------------------------------------
