CFD Online Discussion Forums - Flow in a Diverging Nozzle

Hi
I'm Writing a code in matlab to solve the euler equation for flow in a diverging nozzle with supersonic inlet and subsonic outlet (sectioan 12-5 of Computatioanl fluid dynamics Vol.2 Hoffman). I'm using Steger and warming Flux Vector Splitting explicit formulation. the code doesn't work. Here is the code:
% Allocating variables:
dx = 0.1;
dt = 0.00001;
x = 0:dx:1;
s = 1.398 + 0.374 * tanh(0.8*x - 4);
ds = 0.2776 * ( 1 - tanh(0.8*x - 4).^2 );
gama = 1.4;
n = size(x,2);
ro = zeros(1,n);
u = zeros(1,n);
P = zeros(1,n);
et = zeros(1,n);
Q = zeros(1,n);
H = zeros(1,n);
E_plus = zeros(3,n);
E_minus = zeros(3,n);
DELP = 1;
Iteration = 0;
% Setting initial and boundary conditions for Supersonic inlet and
% subsonic outlet:
u(1:29) = 1676.55;
u(30:n) = 572.76;
ro( : ) = 0.002241;
P( : ) = 2000;
et= 0.5 * u.^2 + P ./ ( ro .* (gama-1) );
a = sqrt(gama.*P./ro);
Q(1,: ) = ro;
Q(2,: ) = ro .* u;
Q(3,: )= ro .* et;
% Sloving the euler equation using Steger and Warming Flux Vector
% Splitting:
while DELP > 0.1
P0 = P;
H(2,: ) = ds .* P;
E_plus(1,: ) = s .* ro .* ( 2 * gama * u + a - u ) / 2 / gama;
E_plus(2,: ) = s .* ro .* ( ( 2 * ( gama - 1 ) * u.^2 ) + ( u + a ).^2 ) / 2 / gama;
E_plus(3,: ) = s .* ro .* ( ( ( gama - 1 ) * u.^3 ) + ( ( u + a ).^3 / 2 )...
+ ( ( ( 3 - gama ) .* ( u + a ) .* a.^2 ) / ( 2 * (gama - 1) ) )) / 2 / gama;
E_minus(1,: ) = s .* ro .* ( u - a ) / 2 / gama;
E_minus(2,: ) = s .* ro .* ( ( u - a ).^2 ) / 2 / gama;
E_minus(3,: ) = s .* ro .* ( ( ( ( u - a ).^3) / 2 ) + ...
( ( 3 - gama ) * ( u - a ) .* a.^2 ) / ( 2 * ( gama - 1 ) ) ) / 2 / gama;
for ii = 2:n-1
if u(ii) > a(ii)
E_plus(:,ii-1:ii) = E_plus(:,ii-1:ii) + E_minus(:,ii-1:ii);
E_minus(:,ii:ii+1) = 0;
end
Q(:,ii) = Q(:,ii) + dt * ( H(:,ii) - ( E_plus(:,ii) - E_plus(:,ii-1)...
+ E_minus(:,ii+1) - E_minus(:,ii) ) / dx ) / s(ii);
end
Q(1,n) = Q(1,n-1);
Q(3,n) = Q(3,n-1);
Q(2,n) = u(n) * ro(n);
ro = Q(1,: );
u( : ) = Q(2,: ) ./ Q(1,: );
et( : ) = Q(3,: ) ./ Q(1,: );
P = ro .* ( gama - 1 ) .* ( et - 0.5 * u.^2);
a = sqrt(gama * P./ro);
DELP = sum( abs(P - P0) );
Iteration = Iteration + 1;
end

What is wrong?
thanks.