How to introduce 2nd order derivative term in pdepe?

Hi,

Refer to the subject cited above, I'm not sure if it is possible or not to introduce second order derivative term in pdepe?

As shown in the image as above, I tried to introduce this 2nd derivative term by a new variable u4 in the code. Such that

u4 = DuDx(3), with source as DuDx(3) and zero flux. The ic and bc are also introduced in the code.

But, the code doesn't work as it gives "Spatial discretization error".

% solve 3-F bifurcaiton model 
function pde2fshear_v4Perturbed_nonlinear
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 sigma_turb;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global z;


  
%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.1 ;
chi_ano=10;
D_ano=10;



% For nonlinear model
gamma_nonlin = 1;
alpha_nonlin =1;
D_0 = 1;
group_vel=0;




%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);


%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 = gradient of turbulence intensity

u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3); 
u4 = sol(:,:,4);

%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;


 % --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global xstep;
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 z;


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



c = [1;1;1;0];
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)=DuDx(3); 


term1 = abs(DuDx(1))-data.constant.critgradpressure;
drift_velFluct = D_0*DuDx(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
        
 

%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*DuDx(3)+D_0*(1+alpha_nonlin)*(alpha_nonlin*u(3)^(alpha_nonlin-1)*(DuDx(3))^2+u(3)^alpha_nonlin*u(4))) ;

 
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2;s3;DuDx(3)];


f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ;D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin;0].*DuDx;                 % flux term for nonlinear model





% --------------------------------------------------------------
function u0 = pdex2ic(x)
 %u0 = [eps; eps; eps;];
 %u0 = [0.01; 0.01; 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);0.5*exp(-100*(x-1)^2)];
% --------------------------------------------------------------
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); 0; ur(4)];
qr = [0.01; 0.1; 1; 1];


%---------------------------------------------------------------

It would help me if someone advice me.

with rgds,

rc

Answer 1

[Torsten]

If all the second-order spatial derivatives in your equation can be written in one expression as d/dx(f(x,t,I,dI/dx)*dI/dx) with a suitably chosen function f, "pdepe" is able to handle your equation. If not, "pdepe" is not able to handle your equation.

Even the term d^2/dx^2(D_I(I)*I) cannot be used without further manipulations within "pdepe" because it's not written in the necessary format.

[Rahul]

Hi Torsten,

Thanks a lot. I applied your suggestion and it worked by combining it with the flux term. The code runs but still one issue remains.

The issue is that although the code runs but it does not converge. I tried it different combinations of input values, but I'm not able to get convergence.

May you please suggest.

with rgds,

rc

[Torsten]

If you include your rewritten equations and boundary conditions and your rewritten code, I have at least a chance to suggest something.

[Rahul]

Hi Torsten,

Please find the code as below which doesn't converge for any set of input values. Also my equation for u3 now looks after solving and rearranging for flux term as below:

% solve 3-F bifurcaiton model 
function pde2fshear_v4Perturbed_nonlinear
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;



  
%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.5;
chi_ano=10;
D_ano=10;



% For nonlinear model
gamma_nonlin = 1;
alpha_nonlin =1;
D_0 = 1;
group_vel=0.5;




%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);


%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));
                                      
        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));
                                   
        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));                       
        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));                       
        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));                        
        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));                       
        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));

        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));
                                    
        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));                     
        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));                        
        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));                       
        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));                       
        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;

    
 % --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global xstep;
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;



%lf FFf F  Fb vfDdength = 0.01 or 1
length=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;


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
        
 

%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+D_0*(1+alpha_nonlin)*alpha_nonlin*u(3)^(alpha_nonlin-1)*DuDx(3))*DuDx(3);

 
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

 



% --------------------------------------------------------------
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)];
% --------------------------------------------------------------
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];


%---------------------------------------------------------------

[Torsten]

Also my equation for u3 now looks after solving and rearranging for flux term as below.

I cannot identify "f" in your expression for dI/dt.

If you think that "f" is the factor in front of d^2I/dx^2, you are mistaken because the factor depends on I and x. Thus d/dx (factor*dI/dx) does not equal factor*d^2I/dx^2.

[Rahul]

Hi Torsten,

I did it as below,

which I coded as in the previous reply of mine.

with rgds,

rc

[Torsten]

Didn't you read my response ?

d/dx((D_I(I) - D_0*(1+alpha)*I^alpha) dI/dx) 

does not equal

(D_I(I) - D_0*(1+alpha)*I^alpha)*d^2I/dx^2

because

D_I(I)-D0*(1+alpha)*I^alpha

depends on x.

Thus

d/dx(D_I(I) - D_0*(1+alpha)*I^alpha dI/dx) = 
(D_I(I) - D_0*(1+alpha)*I^alpha)*d^2I/dx^2 + d/dx(D_I(I) - D_0*(1+alpha)*I^alpha) * dI/dx

[Rahul]

Hi Torsten,

Thanks a lot. Yes, I read your previous response, but didn't implement effectively. Doing now.

with rgds,

rc

[Rahul]

Hi Torsten,

I tried your suggestion. But I got this error

Error using pdepe

Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

The code is as attached.

The equation used is

In the above expression, the portion highlighted in red is implemented using the new variable u4 and the portion in green along with the negative sign is considered as the flux.

plz suggest.

% solve 3-F bifurcaiton model 
function pde2fshear_v4Perturbed_nonlinear
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;



  
%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.5;
chi_ano=10;
D_ano=10;



% For nonlinear model
gamma_nonlin = 1;
alpha_nonlin =1;
D_0 = 5;
group_vel=1;




%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);


%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 =
% D_I-D_0*(1+alpha_nonlin)*I^alpha_nonlin*DuDx(3)


u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3); 
u4 = sol(:,:,4);

%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

    
 % --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global xstep;
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;



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



c = [1;1;1;0];
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;
drift_velFluct = D_0*DuDx(3)^alpha_nonlin;
intensity_diff = D_0*u(3)^alpha_nonlin;
u(4)=intensity_diff-(D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin)*DuDx(3);

% 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 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)+DuDx(4); 

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


f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ;-(intensity_diff-D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin);-D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin].*DuDx;                 % flux term for nonlinear model

 



% --------------------------------------------------------------
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);0.1*(1-x^2)];
% --------------------------------------------------------------
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);0; 0];
qr = [0.01; 0.1; 1; 1];


%---------------------------------------------------------------

[Torsten]

I tried your suggestion. But I got this error

Sorry, but what I wanted to say with my remarks was: your problem is simply not solvable with "pdepe".