code

macnozzle2.m

%NAMED CONSTANTS
RHO = 1;
T = 2;
V = 3;
var = 3;
n_points = 31;
n_time = 14000;

sol = zeros(var, n_points, n_time);
dx = .1;
%initiliaze at t = 0
% A = 1 + 2.2(x-1.5)^2 

C = .02;
X = (0:n_points -1)./10;
A = 1+2.2.*(X-1.5).^2;

for i = 1:n_points
    sol(RHO,i,1) = 1 - .3146.*X(i); 
    sol(T,i,1) = 1 - .2314.*X(i);
    sol(V,i,1) = (.1 + 1.09.*X(i)).*(sol(T,i,1).^.5);
    
end

dts = zeros(1,n_points);

drho1 = zeros(1,n_points);
dV1 = zeros(1,n_points);
dT1 = zeros(1,n_points);
drho2 = zeros(1,n_points);
dV2 = zeros(1,n_points);
dT2 = zeros(1,n_points);
gamma = 1.4;


for t = 1:n_time-1
    %predictor step
    for i = 2:n_points-1 
        
        drho1(i) = -sol(RHO,i,t).*(sol(V,i+1,t)-sol(V,i,t))./dx - sol(RHO,i,t).*sol(V,i,t).*(log(A(i+1)) - log(A(i)))./dx - sol(V,i,t).*(sol(RHO,i+1,t)-sol(RHO,i,t))./dx;
        dV1(i) = -sol(V,i,t).*(sol(V,i+1,t)-sol(V,i,t))./dx - (1./gamma).*( (sol(T,i+1,t) - sol(T,i,t))./dx + (sol(T,i,t)./sol(RHO,i,t)).*(sol(RHO,i+1,t)-sol(RHO,i,t))./dx     );
        dT1(i) = -sol(V,i,t).*(sol(T,i+1,t)-sol(T,i,t))./dx - (gamma -1).*(sol(T,i,t)).*( (sol(V,i+1,t) - sol(V,i,t))./dx +sol(V,i,t).*(log(A(i+1)) - log(A(i)))./dx   );
        
        dts(1,i) = C.*(dx./(sol(T,i,t).^.5 + sol(V,i,t)));
    end
    dt = min(dts(1,2:end-1));
    
    for i = 2:n_points-1
        sol(RHO,i,t+1) = sol(RHO,i,t) + drho1(i).*dt;
        sol(V,i,t+1) = sol(V,i,t) + dV1(i).*dt;
        sol(T,i,t+1) = sol(T,i,t) + dT1(i).*dt;
    end
   
    sol(V,1,t+1)  = sol(V,1,t);
    sol(RHO,1,t+1) = sol(RHO,1,t);
    sol(T,1,t+1) = sol(T,1,t);
    
    sol(V,n_points,t+1)  = sol(V,n_points,t);
    sol(RHO,n_points,t+1) = sol(RHO,n_points,t);
    sol(T,n_points,t+1) = sol(T,n_points,t);
    
    
    
    
    
    
    for i = 2:n_points-1
        drho2(i) = -sol(RHO,i,t+1).*(sol(V,i,t+1)-sol(V,i-1,t+1))./dx - sol(RHO,i,t+1).*sol(V,i,t+1).*(log(A(i)) - log(A(i-1)))./dx - sol(V,i,t+1).*(sol(RHO,i,t+1)-sol(RHO,i-1,t+1))./dx;
        dV2(i) = -sol(V,i,t+1).*(sol(V,i,t+1)-sol(V,i-1,t+1))./dx - (1./gamma).*( (sol(T,i,t+1) - sol(T,i-1,t+1))./dx + (sol(T,i,t+1)./sol(RHO,i,t+1)).*(sol(RHO,i,t+1)-sol(RHO,i-1,t+1))./dx     );
        dT2(i) = -sol(V,i,t+1).*(sol(T,i,t+1)-sol(T,i-1,t+1))./dx - (gamma -1).*(sol(T,i,t+1)).*( (sol(V,i,t+1) - sol(V,i-1,t+1))./dx +sol(V,i,t+1).*(log(A(i)) - log(A(i-1)))./dx   );
        
        drho = mean([drho1(i) drho2(i)]);
        dV = mean([dV1(i) dV2(i)]);
        dT = mean([dT1(i) dT2(i)]);
        
        sol(RHO,i,t+1) = sol(RHO,i,t) + drho.*dt;
        sol(V,i,t+1) = sol(V,i,t) + dV.*dt;
        sol(T,i,t+1) = sol(T,i,t) + dT.*dt;
        
    end
    
    sol(V,1,t+1)  = 2.*sol(V,2,t+1) - sol(V,3,t+1);
    sol(RHO,1,t+1) = sol(RHO,1,t);
    sol(T,1,t+1) = sol(T,1,t);
    
    
    sol(V,n_points,t+1)  = 2.*sol(V,n_points-1,t+1) - sol(V,n_points-2,t+1);
    sol(RHO,n_points,t+1)  = 2.*sol(RHO,n_points-1,t+1) - sol(RHO,n_points-2,t+1);
    sol(T,n_points,t+1)  = 2.*sol(T,n_points-1,t+1) - sol(T,n_points-2,t+1);
    
    
   
    fprintf(num2str(t) + "\n");
end

change the time step(line 46) to dt = min( c*dx ./ (sqrt( sol(T, : ,t) ) + sol(V,:,t) ) )
replace the line 46 to line 39