Why do I have an error in pdepe solving?

Hi all,

I'm having an error while solving a pdepe. It gives me the error

"Index exceeds the number of array elements. Index must not exceed 4."

It will help me immensely if someone kindly help me in understanding as why is it here.

I know in general it happens if array exceeds the no.

But precisely, this seems not the case here in mine.

kindly find the code as attached.

pde2fshear_v4Perturbed_randomWalk()

drift_vel = 
   -0.1000

dCoeff = 
    0.0300

Index exceeds the number of array elements. Index must not exceed 4.

Error in solution>pdex2pde (line 544)
s(4) = chi_growth*(term1*theta_heaviside1-gamma_nonlin*u(5)^alpha_nonlin-alpha_0*u(3)-alpha_v*u(4))*u(5) -drift_vel*DuDx(5);

Error in pdepe (line 242)
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

Error in solution>pde2fshear_v4Perturbed_randomWalk (line 152)
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

% ****************************************************************
% *  Program: pde2fshear_v4Perturbed (Random walk model)           *
% *  Description: Solving coupled thermal, particle transport &  *
% *  turbulence intensity equations, suppressed by flow shear 
% *  and zonal shear                                             *
% *  Author: Rahul Chakraborty                                   *
% *  Email:  [e-mail]                          *
% *  Date: 2023-04-29                                       	 *
% *  Version: 1.0                                          		 *
% ****************************************************************
% solve 3-F bifurcation model 
function pde2fshear_v4Perturbed_randomWalk
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 theta_heaviside1;
global chi_growth;
global beta_growth;        % strength of non-linear growth
global delta_suppress;
global lambda_suppress;
global v_e;
global chi_ano;
global D_ano;
global sigma_turb;



%constants for random walk model
global gamma_nonlin;
global alpha_nonlin;
global drift_vel;
global dCoeff;
global rho_0;
global gamma_damp;
global alpha_0;
global alpha_v;


%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 = 25;  S0 = 21;
track = 1;
chi_growth=0.5;
lambda_suppress=0.01;
sigma_turb=1.5;



gamma_nonlin = 10;
alpha_nonlin = 1;
rho_0=1;
gamma_damp=5;
alpha_0=5;
alpha_v=1;


%position and time grids information
xstep = 100;
tstep = 100; 
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 10;
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);
grad_u4 = zeros(tstep,xstep);
grad_u5 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
curve_u4 = zeros(tstep,xstep);
curve_u5 = 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);


%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;
 %Turbulence intensity equation for Random Walk

  numberofWalk = 5;
  stepsperWalk = 5;
  std_step = 0.1;
  mean_wait = 1.0;                                   %set as appropriate
  delta_X = sign(randn(1,stepsperWalk)) * std_step;
  delta_T = sign(rand(1,stepsperWalk)) * mean_wait;
  
  t_total = transpose(cumsum([0, delta_T]));
  x_total = transpose(cumsum([0, delta_X]));
  
  stem(t_total, x_total)
  xlabel('time'); ylabel('displacement');

  xt = zeros(stepsperWalk,2);
   
  for walk = 1:numberofWalk  
        for step = 1:stepsperWalk
            xt(step+1,1) = xt(step,1) + x_total(step);
            xt(step+1,2) = xt(step,2) + t_total(step);
        end
  end  
 
 plot(xt(:,2),xt(:,1))
 drift_vel = mean(xt(step,1)/xt(step,2));
 dCoeff = mean(xt(step,1)^2/(2*xt(step,2)));
 disp('drift_vel = '); disp(drift_vel);
 
 disp('dCoeff = '); disp(dCoeff);
 
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 = zonal flow
% fourth                   as u4 = mean flow
% fifth                    as u5 = turb intensity
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3); 
u4 = sol(:,:,4);
u5 = sol(:,:,5);

%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));
            grad_u5(j,i) = (u5(j,2)-u5(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));
            grad_u5(j,i) = (u5(j,i)-u5(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(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));
            grad_u5(j,i) = (u5(j,i)-u5(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(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;
         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)) + chi1*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j))*(-grad_u1(i,j)); 
            Q0(i,j) = -grad_u1(i,j)*chi0;
             Q1(i,j) = (chi1*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(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) = chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u5(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;
         else
             Q(i,j) = (chi0*(-grad_u1(i,j)) + chi1*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j))*(-grad_u1(i,j)); 
             Q0(i,j) = -grad_u1(i,j)*chi0;
             Q1(i,j) = (chi1*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(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) = chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
             Gam(i,j) = -grad_u2(i,j)*(D0 + D1*(-grad_u2(i,j)+data.constant.critgraddensity)*u5(i,j)/flowshear_n(i,j));
             Gam0(i,j) = -grad_u2(i,j)*D0;
             Gam1(i,j) = -grad_u2(i,j)*D1*(-grad_u2(i,j)+data.constant.critgraddensity)*u5(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) = D1*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(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));
            curve_u5(j,i) = (grad_u5(j,2)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(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.zonalflow = u3;
data.variable.meanflow = u4;
data.variable.turbulenceintensity = u5;
data.variable.gradpressure = grad_u1;
data.variable.graddensity = grad_u2;
data.variable.gradzonalflow = grad_u3;
data.variable.gradmeanflow = grad_u4;
data.variable.gradturbulenceintensity = grad_u5;
data.variable.curvepressure = curve_u1;
data.variable.curvedensity = curve_u2;
data.variable.curvezonalflow = curve_u3;
data.variable.curvemeanflow = curve_u4;
data.variable.curveturbulenceintensity = curve_u5;
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.predicted.Q = Q_pre;
%data.predicted.Gam = Gam_pre;
%data.predicted.gradientp = pp_pre;
%data.predicted.gradientn = nn_pre;
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 chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global track;
global xstep;
global alpha_chi;
global alpha_D;
global chi_growth;         % total growth rate
global beta_growth;        % strength of non-linear growth
global delta_suppress;
global length;
global tau_conf;
global I0;      % edge localized source of turbulence intensity    
global theta_heaviside1;


%global count;
global v_e;
global drift_vel;


%constants for random walk model
global gamma_nonlin;
global alpha_nonlin;
global drift_vel;
global dCoeff;
global alpha_0;
global rho_0;
global gamma_damp;
global alpha_v;


c = [1;1;1;1];
length=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;

u(4) = mean(v_e)^2;

term1 = abs(DuDx(1))-data.constant.critgradpressure;


% Implementing Heaviside function for H-mode in p, n and I equations
       
if term1 > 0
    theta_heaviside1=1;
else
    theta_heaviside1=0;
end
        
 
 s(3) = -gamma_damp*u(3);
 
 
 s(4) = chi_growth*(term1*theta_heaviside1-gamma_nonlin*u(5)^alpha_nonlin-alpha_0*u(3)-alpha_v*u(4))*u(5) -drift_vel*DuDx(5);
 
 
 s = [H0*exp(-100*x^2/length)+H0/2; S0*exp(-100*(x-0.9)^2/length)+S0/2; s(3);s(4)];

 f = [chi0+chi1*u(5)/flowshear_p; D0+D1*u(5)/flowshear_n; alpha_0*u(3)*u(5)/(1+rho_0*u(4)); dCoeff].*DuDx;                          % flux term for CTRW       


end
% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01; 0.1*exp(-100*(x-1)^2)];
 u0 = [0.4*(1-x^2); 0.4*(1-x^2);0.4*(1-x^2);0.1*exp(-100*(x-1)^2)];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)

pl = [0; 0; 0; 0];
ql = [1; 1; 1; 1];


pr = [ur(1); ur(2); ur(3); 0];
qr = [0.01; 0.001; 1; 1;];

% qr = [0; 0.1; 1];
%qr = [0; eps; 1];
end
%---------------------------------------------------------------

with regards,

rc

Answer 1

[Shivam]

Hi Rahul,

From the information provided, I understand that while executing the attached code, you are getting "index out of bound" issue.

You can debug in your code to know that the error is in line 544 because of indexing in vector 'u' for index 5 whereas the size of vector 'u' is 4x1 only.

>> u

u =
    0.4000
    0.4000
    0.4000
    0.0000

I hope this helps.

Thanks,

Shivam

[Rahul]

Hi,

Thanks.

But the error is only in case of u(5). The program doesn't accept u(5).

I'm unable to figure it out.

rgds,

rc

[Torsten]

You solve for 4 unknown functions. Thus u(5) does not exist.

[Rahul]

Hi Torsten,

Thanks for your comment.

I did changes in the code. But still the error continues.

Attached herewith the code.Can u plz comment.

pde2fshear_v4Perturbed_randomWalk()

Warning: Failure at t=1.855077e-05. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.421011e-20) at time t.

Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00.

Index in position 1 exceeds array bounds. Index must not exceed 1.

Error in solution>pde2fshear_v4Perturbed_randomWalk (line 135)
grad_u1(j,i) = (u1(j,2)-u1(j,1))/(x(2)-x(1));

function pde2fshear_v4Perturbed_randomWalk
global chi0;   % declare global variables
global D0;
global chi1;
global D1;
global alpha_chi;
global alpha_D;
global H0;
global S0;
global Y0;
global H;
global S;
global data;
global track;
global xstep;
global tstep;
global count;
global theta_heaviside1;
global chi_growth;
global lambda_suppress;
global v_e;
global chi_ano;
global D_ano;
global sigma_turb;
global gamma_nonlin;
global alpha_zero;
global alpha_v;
global rho_0;
global gamma_damp;
global alpha_nonlin;
global dCoeff;


  
%define values for constants
data.constant.critgradpressure =  1.2; 
data.constant.critgraddensity =   1.0;
count = 1;

chi0 = 0.5;
D0 = 0.5;
chi1 = 0.5;
D1 = 0.5;
alpha_chi = 0.1;
alpha_D = 0.1;
H0 = 4;
S0 = 13;
track = 1;
chi_growth=20;
lambda_suppress=0.5 ;
sigma_turb=2.0;
chi_ano=10;
D_ano=10;




gamma_nonlin = 10;
gamma_damp = 1;
alpha_zero = 1;
alpha_v = 1;
rho_0 = 1;
alpha_nonlin = 1;
dCoeff = 1;



%position and time grids information
xstep = 100;
tstep = 100; 
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 10;
%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);
grad_u4 = zeros(tstep,xstep);
grad_u5 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
curve_u4 = zeros(tstep,xstep);
curve_u5 = 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);


%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
% fourth solution component as u4 = zonal flow
% fifth solution component as u5 = mean flow
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3); 
u4 = sol(:,:,4);
u5 = sol(:,:,5);



%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));
            grad_u5(j,i) = (u5(j,2)-u5(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));
            grad_u5(j,i) = (u5(j,i)-u5(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(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));
            grad_u5(j,i) = (u5(j,i)-u5(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(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;
         elseif 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+chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(i,j));
            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) = chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(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;
        % elseif 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+D1*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(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);
         else
             % Q(i,j) = -grad_u1(i,j)*(chi0+chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(i,j));
             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) = -grad_u1(i,j)*chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(i,j);
             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) = chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(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+D1*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(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);
         end
     end
 end


%for i=1:xstep
 %   a5(i) = sol_dirac(i) ; 
  %  data.variable.a5_term(i) = a5(i);
%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));
            curve_u5(j,i) = (grad_u5(j,2)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i)-grad_u5(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(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.zonalflow = u4;
data.variable.meanflow = u5;
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.predicted.Q = Q_pre;
%data.predicted.Gam = Gam_pre;
%data.predicted.gradientp = pp_pre;
%data.predicted.gradientn = nn_pre;
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 chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global track;
global xstep;
global tstep;
global alpha_chi;
global alpha_D;
global chi_growth;         % total growth rate
global length;
global theta_heaviside1;
global v_e;
global chi_ano;
global D_ano;
global sigma_turb;
global gamma_nonlin;
global alpha_zero;
global alpha_v;
global rho_0;
global gamma_damp;
global alpha_nonlin;
global dCoeff;

%lf FFf F  Fb vfDdength = 0.01 or 1
length=1;



c = [1;1;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;



term1 = abs(DuDx(1))-data.constant.critgradpressure;



% Implementing Heaviside function for H-mode in p, n and I equations
       
if term1 > 0
    theta_heaviside1=1;
else
    theta_heaviside1=0;
end
        
 

%Turbulence intensity Equation for random walk

s_1 = alpha_zero*u(4)*u(3)/(1+rho_0*u(5))-gamma_damp*u(4);

s_2 = mean(v_e)^2;

s_3 = chi_growth*(term1*theta_heaviside1-gamma_nonlin*u(3)^alpha_nonlin-alpha_zero*u(4)-alpha_v*u(5))*u(3) ;

 s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2; s_1;s_2;s_3];





f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ; 0; 0; dCoeff].*DuDx;                 % flux term for random walk model



%data.function.c(track) = c;
%data.function.f(track) = f;

if track <= xstep
    data.function.heatsource(track) = s(1);
    data.function.particlesource(track) = s(2);
    data.function.turbulencesource(track)= s(3);
    data.function.x(track) = x;
    track = track+1;
end
end

% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01; 0.1*exp(-100*(x-1)^2)];
u0 = [0.4*(1-x^2); 0.4*(1-x^2); 0.5*exp(-100*(x-1)^2);0.4*(1-x^2);0.4*(1-x^2)];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)

pl = [0; 0; 0; 0; 0];
ql = [1; 1; 1; 1; 1];


pr = [ur(1); ur(2); ur(3); ur(4);0];
qr = [0.01; 0.1; 1;1;1];

% qr = [0; 0.1; 1];
%qr = [0; eps; 1];
%---------------------------------------------------------------
end

[Torsten]

Your code doesn't produce syntax errors up to the point where you call "pdepe". The solver starts integrating, but stops after 1.855e-5 s because it cannot continue. This is usually an indication that either the equations are hard to solve or you made some errors in the setup.

I suggest you include the mathematical formulation of the equations you are trying to solve with "pdepe".