I couldnt find where there is error in my logic,although my code is running,I compared an available solution for it with my code,but still couldnt find where logic is wrong?

%From JD Anderson textbook for Computational Fluid Dynamics Chap 8,using
%Maccormack technique of space marching PM wave expansion is solved
%numerically
x=0;
spacestep=0;
y=linspace(0,1,41);%This is the computational y domain

Co=0.5;%Courant no
prompt=input('Enter the coeffcient for artificial viscosity:');
Cy=prompt;
digits(4);
for i=1:41 %Boundary Conditions for x=0 is defined here
    M(i)=2;
    p(i)=101325;
    T(i)=286;
    D(i)=1.225;
    a=(1.4*287*T(i))^0.5;
    u(i)=a*2;
    v(i)=0;
    F1(i)=D(i)*u(i);
    F2(i)=D(i)*(u(i)^2)+p(i);
    F3(i)=F1(i)*v(i);
    F4(i)=3.5*p(i)*u(i)+(F1(i)*(u(i)^2)/2);
    G1(i)=0;
    G2(i)=0;
    G3(i)=p(i);
    G4(i)=0;
end
net=[transpose(D) transpose(u) transpose(v) transpose(p) transpose(T) transpose(M)];
disp(net);
while(spacestep<50)
  if x>10 
     h=40+(x-10)*tan(deg2rad(5.352));%Height of the domain which is a function of x
     def=5.352;%Deflection Angle
  else
     h=40;
     def=0;
  end
  %Loop for calculating the marching step size of the domain across x
  %direction
  for i=1:41
    MA=rad2deg(asin(1/M(i)));
    pos=abs(tan(deg2rad(5.352+MA)));
    neg=abs(tan(deg2rad(5.352-MA)));
    a1=max(pos,neg);
    delX(i)=Co*h*0.025/a1;
  end
  
  delX;
  step=min(delX);%Finding the minimum value of step size of the 41 grid pts along y axis
  [dF1_dx,dF2_dx,dF3_dx,dF4_dx,F1_p,F2_p,F3_p,F4_p,p_p,G1_p,G2_p,G3_p,G4_p]=predictor(F1,F2,F3,F4,G1,G2,G3,G4,def,h,Cy,step,p,y);
  [dF1c_dx,dF2c_dx,dF3c_dx,dF4c_dx,SF1c,SF2c,SF3c,SF4c]=corrector(F1_p,F2_p,F3_p,F4_p,p_p,G1_p,G2_p,G3_p,G4_p,h,Cy,y,def);
  dF1avg_dx=(dF1_dx+dF1c_dx)/2;
  dF2avg_dx=(dF2_dx+dF2c_dx)/2;
  dF3avg_dx=(dF3_dx+dF3c_dx)/2;
  dF4avg_dx=(dF4_dx+dF4c_dx)/2;
  for i=1:41
      F1(i)=F1(i)+SF1c(i)+dF1avg_dx(i)*step;%F1,F2,F3,F4 are the flux terms required to be found
      F2(i)=F2(i)+SF2c(i)+dF2avg_dx(i)*step;
      F3(i)=F3(i)+SF3c(i)+dF3avg_dx(i)*step;
      F4(i)=F4(i)+SF4c(i)+dF4avg_dx(i)*step;
      A=((F3(i)^2)/(2*F1(i)))-F4(i);
      B=3.5*F1(i)*F2(i);
      C=-3*(F1(i)^3);
      D(i)=(-B+(B^2-4*A*C)^0.5)/(2*A);%Here the decoding of primitive variables are done
      u(i)=F1(i)/D(i);
      v(i)=F3(i)/F1(i);

      p(i)=F2(i)-F1(i)*u(i);
      T(i)=p(i)/(D(i)*287);
      magV=(u(i)^2+v(i)^2)^0.5;
      M(i)=magV/((1.4*287*T(i))^0.5);
      G1(i)=D(i)*v(i);%G1,G2,G3,G4 are the source terms
      G2(i)=F3(i);
      G3(i)=D(i)*(v(i)^2)+p(i);  
      G4(i)=3.5*p(i)*v(i)+F1(i)*(v(i)^2)/2;
  end
  v
  
  %Here the Abbot Boundary Condition has been employed for the wall grid point 
  if x>10
       ang=rad2deg(atan(abs(v(1))/u(1)));   
       Rot_ang=def-ang;
  else
       ang=rad2deg(atan(v(1)/u(1)));  
       Rot_ang=ang;
  end
  M
  a=M(1);
  pmf=prandtlmeyer(a);
  pmactual=pmf+Rot_ang;
  act=Mach(1.4,pmactual);
  act_p=p(1)*((1+0.2*a^2)/(1+0.2*act^2))^3.5;
  act_T=T(1)*(1+0.2*a^2)/(1+0.2*act^2);
  D(1)=act_p/(287*act_T);
  p(1)=act_p;
  T(1)=act_T;
  v(1)=-u(1)*tan(deg2rad(def));
  M(1)=act;
  %M(1)=((u(1)^2+v(1)^2)^0.5)/((1.4*287*T(1))^0.5);
  F1(1)=D(1)*u(1);
  F2(1)=D(1)*(u(1)^2)+p(1);
  F3(1)=F1(1)*v(1);
  F4(1)=3.5*p(1)*u(1)+F1(1)*(u(1)^2)/2;
  G1(1)=D(1)*v(1);
  G2(1)=F3(1);
  G3(1)=D(1)*(v(1)^2)+p(1);  
  G4(1)=3.5*p(1)*v(1)+F1(1)*(v(1)^2)/2;
  
  M
  
  
    
  
  
  
  %net=[transpose(D) transpose(u) transpose(v) transpose(p) transpose(T) transpose(M)];
  %disp(net);
  x=x+step;
  spacestep=spacestep+1;
  
end
%Predictor function
function [dF1_dx,dF2_dx,dF3_dx,dF4_dx,F1_p,F2_p,F3_p,F4_p,p_p,G1_p,G2_p,G3_p,G4_p] = predictor(F1,F2,F3,F4,G1,G2,G3,G4,def,h,Cy,step,p,y)
  
   
  for i=1:41
     
    deta_dx=(1-y(i))*tan(deg2rad(def))/h;  
    if i~=41  
      dF1_dy=(F1(i)-F1(i+1))/0.025;
      dG1_dy=(G1(i)-G1(i+1))/0.025;
      dF1_dx(i)=deta_dx*dF1_dy+(dG1_dy/h);
      dF2_dy=(F2(i)-F2(i+1))/0.025;
      dG2_dy=(G2(i)-G2(i+1))/0.025;
      dF2_dx(i)=deta_dx*dF2_dy+(dG2_dy/h);
      dF3_dy=(F3(i)-F3(i+1))/0.025;
      dG3_dy=(G3(i)-G3(i+1))/0.025;
      dF3_dx(i)=deta_dx*dF3_dy+(dG3_dy/h);
      dF4_dy=(F4(i)-F4(i+1))/0.025;
      dG4_dy=(G4(i)-G4(i+1))/0.025;
      dF4_dx(i)=deta_dx*dF4_dy+(dG4_dy/h);
      
    else
      dF1_dy=(F1(i-1)-F1(i))/0.025;
      dG1_dy=(G1(i-1)-G1(i))/0.025;
      dF1_dx(i)=deta_dx*dF1_dy+(dG1_dy/h);
      dF2_dy=(F2(i-1)-F2(i))/0.025;
      dG2_dy=(G2(i-1)-G2(i))/0.025;
      dF2_dx(i)=deta_dx*dF2_dy+(dG2_dy/h);
      dF3_dy=(F3(i-1)-F3(i))/0.025;
      dG3_dy=(G3(i-1)-G3(i))/0.025;
      dF3_dx(i)=deta_dx*dF3_dy+(dG3_dy/h);
      dF4_dy=(F4(i-1)-F4(i))/0.025;
      dG4_dy=(G4(i-1)-G4(i))/0.025;
      dF4_dx(i)=deta_dx*dF4_dy+(dG4_dy/h);
    
    end

    if i==1 || i==41
          SF1=0;
          SF2=0;
          SF3=0;
          SF4=0;
    else  %These are the artificial viscosity terms for predictor step
          SF1=(abs(p(i+1)-2*p(i)+p(i-1))*Cy*(F1(i+1)-2*F1(i)+F1(i-1)))/(p(i+1)+2*p(i)+p(i-1));
          SF2=(abs(p(i+1)-2*p(i)+p(i-1))*Cy*(F2(i+1)-2*F2(i)+F2(i-1)))/(p(i+1)+2*p(i)+p(i-1));
          SF3=(abs(p(i+1)-2*p(i)+p(i-1))*Cy*(F3(i+1)-2*F3(i)+F3(i-1)))/(p(i+1)+2*p(i)+p(i-1));
          SF4=(abs(p(i+1)-2*p(i)+p(i-1))*Cy*(F4(i+1)-2*F4(i)+F4(i-1)))/(p(i+1)+2*p(i)+p(i-1));
    end
    F1_p(i)=F1(i)+SF1+dF1_dx(i)*step;
    F2_p(i)=F2(i)+SF2+dF2_dx(i)*step;
    F3_p(i)=F3(i)+SF3+dF3_dx(i)*step;
    F4_p(i)=F4(i)+SF4+dF4_dx(i)*step;
    A=((F3_p(i)^2)/(2*F1_p(i)))-F4_p(i);
    B=3.5*F1_p(i)*F2_p(i);
    C=-3*(F1_p(i)^3);
    D_p(i)=(-B+(B^2-4*A*C)^0.5)/(2*A);
    u_p(i)=F1_p(i)/D_p(i);
   % v_p(i)=F3_p(i)/F1_p(i);
    p_p(i)=F2_p(i)-F1_p(i)*u_p(i);
    %T_p(i)=p_p(i)/(D_p(i)*287);
    G1_p(i)=D_p(i)*F3_p(i)/F1_p(i);
    G2_p(i)=F3_p(i);
    G3_p(i)=D_p(i)*(F3_p(i)/F1_p(i))^2+F2_p(i)-(F1_p(i)^2/D_p(i));
    te=F3_p(i)/F1_p(i);
    G4_p(i)=3.5*(F2_p(i)-(F1_p(i)^2/D_p(i)))*te+0.5*D_p(i)*te*((F1_p(i)/D_p(i))^2+te^2);
    
  end

end
%Corrector function
function [dF1c_dx,dF2c_dx,dF3c_dx,dF4c_dx,SF1c,SF2c,SF3c,SF4c] = corrector(F1_p,F2_p,F3_p,F4_p,p_p,G1_p,G2_p,G3_p,G4_p,h,Cy,y,def)
 
  for i=1:41
    deta_dx=(1-y(i))*tan(deg2rad(def))/h;  
    if i>1  
      dF1_dy=(F1_p(i-1)-F1_p(i))/0.025;
      dG1_dy=(G1_p(i-1)-G1_p(i))/0.025;
      dF1c_dx(i)=deta_dx*dF1_dy+(dG1_dy/h);
      dF2_dy=(F2_p(i-1)-F2_p(i))/0.025;
      dG2_dy=(G2_p(i-1)-G1_p(i))/0.025;
      dF2c_dx(i)=deta_dx*dF2_dy+(dG2_dy/h);
      dF3_dy=(F3_p(i-1)-F3_p(i))/0.025;
      dG3_dy=(G3_p(i-1)-G3_p(i))/0.025;
      dF3c_dx(i)=deta_dx*dF3_dy+(dG3_dy/h);
      dF4_dy=(F4_p(i-1)-F4_p(i))/0.025;
      dG4_dy=(G4_p(i-1)-G4_p(i))/0.025;
      dF4c_dx(i)=deta_dx*dF4_dy+(dG4_dy/h);
    else
      dF1_dy=(F1_p(i)-F1_p(i+1))/0.025;
      dG1_dy=(G1_p(i)-G1_p(i+1))/0.025;
      dF1c_dx(i)=deta_dx*dF1_dy+(dG1_dy/h);
      dF2_dy=(F2_p(i)-F2_p(i+1))/0.025;
      dG2_dy=(G2_p(i)-G1_p(i+1))/0.025;
      dF2c_dx(i)=deta_dx*dF2_dy+(dG2_dy/h);
      dF3_dy=(F3_p(i)-F3_p(i+1))/0.025;
      dG3_dy=(G3_p(i)-G3_p(i+1))/0.025;
      dF3c_dx(i)=deta_dx*dF3_dy+(dG3_dy/h);
      dF4_dy=(F4_p(i)-F4_p(i+1))/0.025;
      dG4_dy=(G4_p(i)-G4_p(i+1))/0.025;
      dF4c_dx(i)=deta_dx*dF4_dy+(dG4_dy/h); 
    
        
    end
    SF1c(1)=0;
    SF2c(1)=0;
    SF3c(1)=0;
    SF4c(1)=0;
    for i=2:40
        SF1c(i)=(abs(p_p(i+1)-2*p_p(i)+p_p(i-1))*Cy*(F1_p(i+1)-2*F1_p(i)+F1_p(i-1)))/(p_p(i+1)+2*p_p(i)+p_p(i-1));
        SF2c(i)=(abs(p_p(i+1)-2*p_p(i)+p_p(i-1))*Cy*(F2_p(i+1)-2*F2_p(i)+F2_p(i-1)))/(p_p(i+1)+2*p_p(i)+p_p(i-1));
        SF3c(i)=(abs(p_p(i+1)-2*p_p(i)+p_p(i-1))*Cy*(F3_p(i+1)-2*F3_p(i)+F3_p(i-1)))/(p_p(i+1)+2*p_p(i)+p_p(i-1));
        SF4c(i)=(abs(p_p(i+1)-2*p_p(i)+p_p(i-1))*Cy*(F4_p(i+1)-2*F4_p(i)+F4_p(i-1)))/(p_p(i+1)+2*p_p(i)+p_p(i-1));
    end
    SF1c(41)=0;
    SF2c(41)=0;
    SF3c(41)=0;
    SF4c(41)=0;
    %dF1c_dx(41)=0;
    %dF2c_dx(41)=0;
    %dF3c_dx(41)=0;
    %dF4c_dx(41)=0;

  end
end
function [pm]=prandtlmeyer(m)
   c=(1.4+1)/(1.4-1);
   
   p2=m^2 -1;
   p3=rad2deg(atan(sqrt(p2/c)));
   p4=rad2deg(atan(sqrt(p2)));
   pm=sqrt(c)*p3-p4;
end
function [t]=Mach(g,p)
 
 b=2.9;
 a=1.1;
 precision=0.0000001; 
 c=(g+1)/(g-1);
 %To find mach,bisection method has been employed
 p2=a^2 -1;
 r2=((a+b)/2)^2-1;
 p3=rad2deg(atan(sqrt(p2/c)));
 r3=rad2deg(atan(sqrt(r2/c)));
 r4=rad2deg(atan(sqrt(r2)));
 p4=rad2deg(atan(sqrt(p2)));
 pm=sqrt(c)*p3-p4;
 rm=sqrt(c)*r3-r4;
 zero1=pm-p;
 zero2=rm-p;
 while((b-a)/2>precision)
      
      if zero1*zero2<=0
          b=(a+b)/2;
      else
          a=(a+b)/2;
      end
      p2=a^2 -1;
      r2=((a+b)/2)^2-1;
      p3=rad2deg(atan(sqrt(p2/c)));
      r3=rad2deg(atan(sqrt(r2/c)));
      r4=rad2deg(atan(sqrt(r2)));
      p4=rad2deg(atan(sqrt(p2)));
      pm=sqrt(c)*p3-p4;
      rm=sqrt(c)*r3-r4;
      zero1=pm-p;
      zero2=rm-p;
      
 end
 t=(a+b)/2;
end


       

[Jan]

In opposite to you, we do not have a valid solution. You do noit mention, what the difference between your result and the example solution is. So all we see is running code, but no hints what it should compute. How could we find a problem?

[Kishore S]

Available solution for Prandtl Meyer Expansion:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NUMERICAL SOLUTION OF A PRANDTL-MEYER EXPANSION WAVE FLOW FIELD %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% By: Iván Padilla ([e-mail]) ; Alberto García ([e-mail])
        
% This program is a numerical solver of the Prandtl-Meyer expansion wave
% problem. It is based on the method provided by Anderson in his book
% "Computational Fluid Dynamics", chapter 8. You may need to take a look at it
% in order to unsderstand the following code.
% Assumptions: inviscid, steady, two-dimensional flow field for a
% calloricaly perfect gas. This simulation is done using standard air at
% Mach 2
close all;
clear all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Geometry and physical domain parameters %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H = 40; % Physical domain height [m]
Theta = 5.352*pi/180; % Expansion angle [rad]
L = 65; % Physical domain longitude [m]
E = 10; % Expansion point distance [m]
%%%%%%%%%%%%%%%%%%%%%
% Inflow conditions %
%%%%%%%%%%%%%%%%%%%%%
M_in = 2; % Inlet Mach number
P_in = 101325; % Inlet pressure (static) [Pa]
Rho_in = 1.225; % Inlet density [kg/m^3]
R_air = 287; % Air gas constant [J/kg*K]
Gamma = 1.4; % Specific heat ratio for perfect air 
T_in = P_in/(Rho_in*R_air); % Inlet temperature [K]
a_in = sqrt(Gamma*R_air*T_in); % Inlet sound speed [m/s]
v_in = 0; % Inlet y-component of velocity [m/s]
u_in = M_in*a_in; % Inlet x-component of velocity [m/s]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flow field variables initialization %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Flow_field = struct('M',zeros(401,877),'u',zeros(401,877),'v',zeros(401,877),'T',zeros(401,877),'P',zeros(401,877),'Rho',zeros(401,877),'a',zeros(401,877),'M_angle',zeros(401,877));
% I have chosen this matrix size in order to preallocate them and increase
% performance. Anderson uses 41 steps in the vertical direction but I
% increased the number of points by a factor of 10 to increase precision.
% This code also works for a grid of 41x87
% Now we are able to fed in the initial data line (first vertical line)
Flow_field.M(:,1) = M_in;
Flow_field.u(:,1) = u_in;
Flow_field.v(:,1) = v_in;
Flow_field.T(:,1) = T_in;
Flow_field.P(:,1) = P_in;
Flow_field.Rho(:,1) = Rho_in;
Flow_field.a(:,1) = a_in;
Flow_field.M_angle(:,1) = asin(1/M_in); % Mach angle [rad]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flux/Divergence terms definition %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = struct('F1',zeros(401,877),'F2',zeros(401,877),'F3',zeros(401,877),'F4',zeros(401,877));
G = struct('G1',zeros(401,877),'G2',zeros(401,877),'G3',zeros(401,877),'G4',zeros(401,877));
F.F1(:,1) = Rho_in*u_in; % From continuity: differentiated with respect to x
G.G1(:,1) = Rho_in*v_in; % From continuity: differentiated with respect to y
F.F2(:,1) = Rho_in*(u_in^2) + P_in; % From x-momentum (Euler): differentiated with respect to x
G.G2(:,1) = Rho_in*u_in*v_in; % From x-momentum (Euler): differentiated with respect to y
F.F3 = G.G2; % From y-momentum (Euler): differentiated with respect to x
G.G3(:,1) = Rho_in*(v_in^2) + P_in; % From y-momentum (Euler): differentiated with respect to y
F.F4(:,1) = (Gamma/(Gamma - 1))*P_in*u_in + Rho_in*u_in*(((u_in^2) + (v_in^2)))/2; % From energy: differentiated with respect to x
G.G4(:,1) = (Gamma/(Gamma - 1))*P_in*v_in + Rho_in*v_in*(((u_in^2) + (v_in^2)))/2; % From energy: differentiated with respect to y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grid transformation and parameters %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A boundary fitted coordinate system is used. The transformation is given
% by the relation between the physical (x,y) and the computational
% (Xi,Eta) planes.
i = 1; % Grid horitzontal position (matrix column)
j = 1; % Grid vertical position (matrix row)
Grid = struct('x',zeros(1,877),'y',zeros(401,877),'y_t',zeros(401,1),'y_s',zeros(1,877),'h',zeros(1,877),'m1',zeros(401,877),'m2',zeros(401,877),'delta_y_t',0.0025);
% x is the horitzontal coordinate position in both the physical and the
% computational planes. So x = Xi
% y is the vertical coordinate position in the physical plane
% y_t is the vertical coordinate position in the computational (transformed) plane (Eta value)
% y_s is the transformation parameter that represents the wall expansion
% geometry
% h is the transformation parameter that represents the vertical distance
% between the wall and the upper boundary
% m1 is the metric representing the derivative of Eta to x in the
% transformation
% m2 is the metric representing the derivative of Eta to y in the
% transformation
for j = 2:401 % Here we compute the vertical component (Eta) of the grid in the computational plane
    Grid.y_t(j) = Grid.y_t(j-1) + Grid.delta_y_t; 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Numerical solution using MacCormack's predictor-corrector technique %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization of variables needed for the technique
% Predicted values (all variables finishing with "_p"
F_p = struct('F1_p',zeros(401,877),'F2_p',zeros(401,877),'F3_p',zeros(401,877),'F4_p',zeros(401,877));
G_p = struct('G1_p',zeros(401,877),'G2_p',zeros(401,877),'G3_p',zeros(401,877),'G4_p',zeros(401,877));
dF_x = struct('dF1_x',zeros(401,1),'dF2_x',zeros(401,1),'dF3_x',zeros(401,1),'dF4_x',zeros(401,1));
dF_p_x = struct('dF1_p_x',zeros(401,1),'dF2_p_x',zeros(401,1),'dF3_p_x',zeros(401,1),'dF4_p_x',zeros(401,1));
P_p = zeros(401,1); % Needed for artificial viscosity in the corrector step
while (Grid.x(i) <= L)
    % First of all compute the grid transformation for the current (i)
    % position:
    Grid = compute_transformation(i,Grid,H,E,Theta);
    
    % Compute the downstream step size (delta_x):
    delta_x = compute_step_size(Theta,i,Grid.y,Flow_field.M_angle);
    
    % Apply predictor step and obtain the derivatives of F with respect to
    % x and the predicted values of F in i+1:
    [dF_x,F_p] = predictor_step(i,F,G,Grid,F_p,dF_x,Flow_field.P,delta_x);
    
    % Obtain the predicted G and P values in i+1:
    [P_p,G_p] = compute_predicted_G(F_p,G_p,P_p,Gamma,i);
    
    % Now we are able to apply corrector step in oder to obtain the desired values
    % of F at the next downstream point i+1:
    F = corrector_step(i,F,F_p,G_p,P_p,dF_x,dF_p_x,Grid,delta_x);
    
    % The last step is to decode the flow field variables at i+1 and
    % compute the new G values. This function also includes the application
    % of Abbett's boundary condition for j = 1:
    [Flow_field,G] = decode_flow_field(i,F,G,Flow_field,Grid,Gamma,R_air,Theta,E);
    
    % Finally advance to the next downstream position
    Grid.x(i+1) = Grid.x(i) + delta_x;
    i = i + 1;
    clc;
    %fprintf('Step %d of 877\n',i);
end
Grid = compute_transformation(i,Grid,H,E,Theta); % To compute the grid transformation at the last point (i = 877)
%disp('Finished!');
mesh(Grid.y,Flow_field.M); % Mesh function plots in 3D, but this has no sense in this simulation. Move the axes to view it in 2d.
xlabel('Downstream position (i)');
zlabel('y (m)');
title('Mach Number');
colorbar;

function [Predicted_pressure_value,Predicted_G_values] = compute_predicted_G(F_p,G_p,P_p,Gamma,i)
for j = 1:401,
    A_p = ((F_p.F3_p(j,i+1)^2)/(2*F_p.F1_p(j,i+1))) - F_p.F4_p(j,i+1);
    B_p = (Gamma/(Gamma - 1))*F_p.F1_p(j,i+1)*F_p.F2_p(j,i+1);
    C_p = -(((Gamma + 1)/(2*(Gamma - 1)))*(F_p.F1_p(j,i+1)^3));
    Rho_p = (-B_p + (sqrt((B_p^2) - (4*A_p*C_p))))/(2*A_p);
    P_p(j) = F_p.F2_p(j,i+1) - ((F_p.F1_p(j,i+1)^2)/Rho_p);
    G_p.G1_p(j,i+1) = Rho_p*(F_p.F3_p(j,i+1)/F_p.F1_p(j,i+1));
    G_p.G2_p(j,i+1) = F_p.F3_p(j,i+1);
    G_p.G3_p(j,i+1) = (Rho_p*((F_p.F3_p(j,i+1)/F_p.F1_p(j,i+1))^2)) + F_p.F2_p(j,i+1) - ((F_p.F1_p(j,i+1)^2)/Rho_p);
    G_p.G4_p(j,i+1) = ((Gamma/(Gamma - 1))*((F_p.F2_p(j,i+1)) - ((F_p.F1_p(j,i+1)^2)/Rho_p))*(F_p.F3_p(j,i+1)/F_p.F1_p(j,i+1))) + (((Rho_p*F_p.F3_p(j,i+1))/(2*F_p.F1_p(j,i+1)))*(((F_p.F1_p(j,i+1)/Rho_p)^2) + ((F_p.F3_p(j,i+1)/F_p.F1_p(j,i+1))^2)));
end
Predicted_pressure_value = P_p;
Predicted_G_values = G_p;
end

function [delta_x] = compute_step_size(Theta,i,y,M_angle)
% Anderson does not give much details about the computation of this step size so
% the following function is how I understood the technique. The obtained
% results are good so I believe this function works well.
tan_max_pos = 0;
tan_max_neg = 0;
for j = 1:401,  % We need to find the maximum value of both tangents, so we need to compute them in the whole vertical line
    tan_pos = abs(tan(Theta+M_angle(j,i)));
    tan_neg = abs(tan(Theta-M_angle(j,i)));
    if (tan_pos > tan_max_pos)
        tan_max_pos = tan_pos;
    end
    if (tan_neg > tan_max_neg)
        tan_max_neg = tan_neg;
    end
end
tan_max = max(tan_max_pos,tan_max_neg);
delta_y = y(2,i) - y(1,i);  % delta_y is constant along the same vertical line (i). 
% The grid transformation geometry imposes this. So it doesnt' matter to
% compute delta_y using y in 2,1 or using it in 21,20 for example, the
% result will be the same.
delta_x = 0.5*(delta_y/tan_max); % Where 0.5 is the Courant number
end

function [Updated_Grid] = compute_transformation(i,Grid,H,E,Theta)
if (Grid.x(i) < E)
    Grid.y_s(i) = 0;
    Grid.h(i) = H;   
else
    Grid.y_s(i) = (-Grid.x(i)*tan(Theta)) + (E*tan(Theta));
    Grid.h(i) = H + (Grid.x(i)*tan(Theta)) - (E*tan(Theta));
end
for j = 1:401
    if (Grid.x(i) < E)
        Grid.m1(j,i) = 0;
    else
        Grid.m1(j,i) = (tan(Theta)/Grid.h(i)) - (Grid.y_t(j)*(tan(Theta)/Grid.h(i)));
    end
    Grid.y(j,i) = (Grid.y_t(j)*Grid.h(i)) + Grid.y_s(i);
    Grid.m2(j,i) = 1/Grid.h(i);
end
Updated_Grid = Grid;
end

function [F_values] = corrector_step(i,F,F_p,G_p,P_p,dF_x,dF_p_x,Grid,delta_x)
for j = 1:401
    if (j == 1) % Forward difference without artificial viscosity (no artificial viscosity in the boundaries)
        dF_p_x.dF1_p_x(j) = (Grid.m1(j,i)*((F_p.F1_p(j,i+1) - F_p.F1_p(j+1,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G1_p(j,i+1) - G_p.G1_p(j+1,i+1))/Grid.delta_y_t)); % The derivative of F1_p with respect to x (Xi) at the point j,i+1
        dF_p_x.dF2_p_x(j) = (Grid.m1(j,i)*((F_p.F2_p(j,i+1) - F_p.F2_p(j+1,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G2_p(j,i+1) - G_p.G2_p(j+1,i+1))/Grid.delta_y_t));
        dF_p_x.dF3_p_x(j) = (Grid.m1(j,i)*((F_p.F3_p(j,i+1) - F_p.F3_p(j+1,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G3_p(j,i+1) - G_p.G3_p(j+1,i+1))/Grid.delta_y_t));
        dF_p_x.dF4_p_x(j) = (Grid.m1(j,i)*((F_p.F4_p(j,i+1) - F_p.F4_p(j+1,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G4_p(j,i+1) - G_p.G4_p(j+1,i+1))/Grid.delta_y_t));
    elseif (j == 401) % Rearward difference without artificial viscosity
        dF_p_x.dF1_p_x(j) = (Grid.m1(j,i)*((F_p.F1_p(j-1,i+1) - F_p.F1_p(j,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G1_p(j-1,i+1) - G_p.G1_p(j,i+1))/Grid.delta_y_t));
        dF_p_x.dF2_p_x(j) = (Grid.m1(j,i)*((F_p.F2_p(j-1,i+1) - F_p.F2_p(j,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G2_p(j-1,i+1) - G_p.G2_p(j,i+1))/Grid.delta_y_t));
        dF_p_x.dF3_p_x(j) = (Grid.m1(j,i)*((F_p.F3_p(j-1,i+1) - F_p.F3_p(j,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G3_p(j-1,i+1) - G_p.G3_p(j,i+1))/Grid.delta_y_t));
        dF_p_x.dF4_p_x(j) = (Grid.m1(j,i)*((F_p.F4_p(j-1,i+1) - F_p.F4_p(j,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G4_p(j-1,i+1) - G_p.G4_p(j,i+1))/Grid.delta_y_t));
    else % Rearward difference with artificial viscosity
        dF_p_x.dF1_p_x(j) = (Grid.m1(j,i)*((F_p.F1_p(j-1,i+1) - F_p.F1_p(j,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G1_p(j-1,i+1) - G_p.G1_p(j,i+1))/Grid.delta_y_t));
        dF_p_x.dF2_p_x(j) = (Grid.m1(j,i)*((F_p.F2_p(j-1,i+1) - F_p.F2_p(j,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G2_p(j-1,i+1) - G_p.G2_p(j,i+1))/Grid.delta_y_t));
        dF_p_x.dF3_p_x(j) = (Grid.m1(j,i)*((F_p.F3_p(j-1,i+1) - F_p.F3_p(j,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G3_p(j-1,i+1) - G_p.G3_p(j,i+1))/Grid.delta_y_t));
        dF_p_x.dF4_p_x(j) = (Grid.m1(j,i)*((F_p.F4_p(j-1,i+1) - F_p.F4_p(j,i+1))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G_p.G4_p(j-1,i+1) - G_p.G4_p(j,i+1))/Grid.delta_y_t));
        
        SF1_p = (((0.6*(abs(P_p(j+1)) - (2*P_p(j)) + P_p(j-1))))/(P_p(j+1) + 2*P_p(j) + P_p(j-1)))*(F_p.F1_p(j+1,i+1) - (2*F_p.F1_p(j,i+1)) + F_p.F1_p(j-1,i+1)); % Predicted value of artificial viscosity at j,i+1
        SF2_p = (((0.6*(abs(P_p(j+1)) - (2*P_p(j)) + P_p(j-1))))/(P_p(j+1) + 2*P_p(j) + P_p(j-1)))*(F_p.F2_p(j+1,i+1) - (2*F_p.F2_p(j,i+1)) + F_p.F2_p(j-1,i+1));
        SF3_p = (((0.6*(abs(P_p(j+1)) - (2*P_p(j)) + P_p(j-1))))/(P_p(j+1) + 2*P_p(j) + P_p(j-1)))*(F_p.F3_p(j+1,i+1) - (2*F_p.F3_p(j,i+1)) + F_p.F3_p(j-1,i+1));
        SF4_p = (((0.6*(abs(P_p(j+1)) - (2*P_p(j)) + P_p(j-1))))/(P_p(j+1) + 2*P_p(j) + P_p(j-1)))*(F_p.F4_p(j+1,i+1) - (2*F_p.F4_p(j,i+1)) + F_p.F4_p(j-1,i+1));
    end
    
    dF1_x_av = 0.5*(dF_x.dF1_x(j) + dF_p_x.dF1_p_x(j));
    dF2_x_av = 0.5*(dF_x.dF2_x(j) + dF_p_x.dF2_p_x(j));
    dF3_x_av = 0.5*(dF_x.dF3_x(j) + dF_p_x.dF3_p_x(j));
    dF4_x_av = 0.5*(dF_x.dF4_x(j) + dF_p_x.dF4_p_x(j));
        
    if (j == 1 || j == 401)
        F.F1(j,i+1) = F.F1(j,i) + (dF1_x_av*delta_x);
        F.F2(j,i+1) = F.F2(j,i) + (dF2_x_av*delta_x);
        F.F3(j,i+1) = F.F3(j,i) + (dF3_x_av*delta_x);
        F.F4(j,i+1) = F.F4(j,i) + (dF4_x_av*delta_x);
    else
        F.F1(j,i+1) = F.F1(j,i) + (dF1_x_av*delta_x) + SF1_p;
        F.F2(j,i+1) = F.F2(j,i) + (dF2_x_av*delta_x) + SF2_p;
        F.F3(j,i+1) = F.F3(j,i) + (dF3_x_av*delta_x) + SF3_p;
        F.F4(j,i+1) = F.F4(j,i) + (dF4_x_av*delta_x) + SF4_p;
    end
end
F_values = F;
end

function [Updated_flow_field,Updated_G] = decode_flow_field(i,F,G,Flow_field,Grid,Gamma,R_air,Theta,E)
for j = 1:401,
    A = ((F.F3(j,i+1)^2)/(2*F.F1(j,i+1))) - F.F4(j,i+1);
    B = (Gamma/(Gamma - 1))*F.F1(j,i+1)*F.F2(j,i+1);
    C = -(((Gamma + 1)/(2*(Gamma - 1)))*(F.F1(j,i+1)^3));
    if (j == 1) % Apply Abbett's wall boundary condition
        Rho_cal = (-B + (sqrt((B^2) - (4*A*C))))/(2*A);
        u_cal = F.F1(j,i+1)/Rho_cal;
        v_cal = F.F3(j,i+1)/F.F1(j,i+1);
        P_cal = F.F2(j,i+1) - (F.F1(j,i+1)*u_cal);
        T_cal = P_cal/(R_air*Rho_cal);
        M_cal = (sqrt((u_cal^2) + (v_cal)^2))/(sqrt(Gamma*R_air*T_cal));
        if (Grid.x(i) < E)
            phi = atan(v_cal/u_cal);
        else
            phi = Theta - (atan(abs(v_cal)/u_cal));
        end
        f_cal = sqrt((Gamma + 1)/(Gamma - 1))*(atan(sqrt(((Gamma - 1)/(Gamma + 1))*(M_cal^2 - 1)))) - (atan(sqrt((M_cal^2) - 1))); % Prandtl-Meyer function
        f_act = f_cal + phi; % Rotated Prandtl-Meyer function
        % We need to find the Mach number corresponding to a Prandtl-Meyer
        % function of value "f_act". Anderson suggests a simple trial and
        % error computation, but I have used a bisection method that I think should be
        % more efficient.
        a_int = 1.1; % Left limit of the interval 
        b_int = 2.9; % Right limit of the interval
        precision = 0.0000001; % Max error 
        zero_f1 = sqrt((Gamma + 1)/(Gamma - 1))*(atan(sqrt(((Gamma - 1)/(Gamma + 1))*(a_int^2 - 1)))) - (atan(sqrt((a_int^2) - 1))) - f_act; % Function used to find its zero
        zero_f2 = sqrt((Gamma + 1)/(Gamma - 1))*(atan(sqrt(((Gamma - 1)/(Gamma + 1))*(((a_int + b_int)/2)^2 - 1)))) - (atan(sqrt((((a_int + b_int)/2)^2) - 1))) - f_act;
        while ((b_int-a_int)/2 > precision)
            if (zero_f1*zero_f2 <=0)
                b_int = (a_int + b_int)/2;
            else
                a_int = (a_int + b_int)/2;
            end
            zero_f1 = sqrt((Gamma + 1)/(Gamma - 1))*(atan(sqrt(((Gamma - 1)/(Gamma + 1))*(a_int^2 - 1)))) - (atan(sqrt((a_int^2) - 1))) - f_act;
            zero_f2 = sqrt((Gamma + 1)/(Gamma - 1))*(atan(sqrt(((Gamma - 1)/(Gamma + 1))*(((a_int + b_int)/2)^2 - 1)))) - (atan(sqrt((((a_int + b_int)/2)^2) - 1))) - f_act;
        end
        M_act = (a_int + b_int)/2; % Corrected Mach number
        Flow_field.M(j,i+1) = M_act;
        Flow_field.M_angle(j,i+1) = (asin(1/Flow_field.M(j,i+1)));
        Flow_field.P(j,i+1) = P_cal*(((1 + ((Gamma - 1)/2)*(M_cal^2))/(1 + ((Gamma - 1)/2)*(M_act^2)))^(Gamma/(Gamma - 1)));
        Flow_field.T(j,i+1) = T_cal*((1 + ((Gamma - 1)/2)*(M_cal^2))/(1 + ((Gamma - 1)/2)*(M_act^2)));
        Flow_field.Rho(j,i+1) = Flow_field.P(j,i+1)/(R_air*Flow_field.T(j,i+1));
        Flow_field.u(j,i+1) = u_cal;
        Flow_field.a(j,i+1) = sqrt(Gamma*R_air*Flow_field.T(j,i+1));
        if (Grid.x(i) > E)
            Flow_field.v(j,i+1) = -(Flow_field.u(j,i+1)*tan(Theta));
        else
            Flow_field.v(j,i+1) = 0;
        end
        % Finally we correct the F terms
        F.F1(j,i+1) = Flow_field.Rho(j,i+1)*Flow_field.u(j,i+1);
        F.F2(j,i+1) = (Flow_field.Rho(j,i+1)*(Flow_field.u(j,i+1)^2)) + Flow_field.P(j,i+1);
        F.F3(j,i+1) = Flow_field.Rho(j,i+1)*Flow_field.u(j,i+1)*Flow_field.v(j,i+1);
        F.F4(j,i+1) = ((Gamma/(Gamma - 1))*Flow_field.P(j,i+1)*Flow_field.u(j,i+1)) + (Flow_field.Rho(j,i+1)*Flow_field.u(j,i+1)*(((Flow_field.u(j,i+1)^2) + (Flow_field.v(j,i+1)^2)))/2);
    else
        Flow_field.Rho(j,i+1) = (-B + (sqrt((B^2) - (4*A*C))))/(2*A);
        Flow_field.u(j,i+1) = F.F1(j,i+1)/Flow_field.Rho(j,i+1);
        Flow_field.v(j,i+1) = F.F3(j,i+1)/F.F1(j,i+1);
        Flow_field.P(j,i+1) = F.F2(j,i+1) - (F.F1(j,i+1)*Flow_field.u(j,i+1));
        Flow_field.T(j,i+1) = Flow_field.P(j,i+1)/(R_air*Flow_field.Rho(j,i+1));
        Flow_field.a(j,i+1) = sqrt(Gamma*R_air*Flow_field.T(j,i+1));
        Flow_field.M(j,i+1) = (sqrt((Flow_field.u(j,i+1)^2) + (Flow_field.v(j,i+1)^2)))/Flow_field.a(j,i+1);
        Flow_field.M_angle(j,i+1) = asin(1/Flow_field.M(j,i+1));
    end
    G.G1(j,i+1) = Flow_field.Rho(j,i+1)*(F.F3(j,i+1)/F.F1(j,i+1));
    G.G2(j,i+1) = F.F3(j,i+1);
    G.G3(j,i+1) = (Flow_field.Rho(j,i+1)*((F.F3(j,i+1)/F.F1(j,i+1))^2)) + F.F2(j,i+1) - ((F.F1(j,i+1)^2)/Flow_field.Rho(j,i+1));
    G.G4(j,i+1) = ((Gamma/(Gamma - 1))*((F.F2(j,i+1)) - ((F.F1(j,i+1)^2)/Flow_field.Rho(j,i+1)))*(F.F3(j,i+1)/F.F1(j,i+1))) + (((Flow_field.Rho(j,i+1)*F.F3(j,i+1))/(2*F.F1(j,i+1)))*(((F.F1(j,i+1)/Flow_field.Rho(j,i+1))^2) + ((F.F3(j,i+1)/F.F1(j,i+1))^2)));
end
Updated_flow_field = Flow_field;
Updated_G = G;
end

function [Updated_derivatives, Predicted_F_values] = predictor_step(i,F,G,Grid,F_p,dF_x,P,delta_x)
for j = 1:401,
    if (j == 1) % Forward difference without artificial viscosity (no artificial viscosity in the boundaries)
        dF_x.dF1_x(j) = (Grid.m1(j,i)*((F.F1(j,i) - F.F1(j+1,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G1(j,i) - G.G1(j+1,i))/Grid.delta_y_t)); % The derivative of F1 with respect to x (Xi) at the point j,i
        dF_x.dF2_x(j) = (Grid.m1(j,i)*((F.F2(j,i) - F.F2(j+1,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G2(j,i) - G.G2(j+1,i))/Grid.delta_y_t));
        dF_x.dF3_x(j) = (Grid.m1(j,i)*((F.F3(j,i) - F.F3(j+1,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G3(j,i) - G.G3(j+1,i))/Grid.delta_y_t));
        dF_x.dF4_x(j) = (Grid.m1(j,i)*((F.F4(j,i) - F.F4(j+1,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G4(j,i) - G.G4(j+1,i))/Grid.delta_y_t));
        F_p.F1_p(j,i+1) = F.F1(j,i) + (dF_x.dF1_x(j)*delta_x); % By the Euler numerical method advance F1 to the next downstream point and obtain its predicted value
        F_p.F2_p(j,i+1) = F.F2(j,i) + (dF_x.dF2_x(j)*delta_x);
        F_p.F3_p(j,i+1) = F.F3(j,i) + (dF_x.dF3_x(j)*delta_x);
        F_p.F4_p(j,i+1) = F.F4(j,i) + (dF_x.dF4_x(j)*delta_x);
    elseif (j == 401) % Rearward difference without artificial viscosity
        dF_x.dF1_x(j) = (Grid.m1(j,i)*((F.F1(j-1,i) - F.F1(j,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G1(j-1,i) - G.G1(j,i))/Grid.delta_y_t));
        dF_x.dF2_x(j) = (Grid.m1(j,i)*((F.F2(j-1,i) - F.F2(j,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G2(j-1,i) - G.G2(j,i))/Grid.delta_y_t));
        dF_x.dF3_x(j) = (Grid.m1(j,i)*((F.F3(j-1,i) - F.F3(j,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G3(j-1,i) - G.G3(j,i))/Grid.delta_y_t));
        dF_x.dF4_x(j) = (Grid.m1(j,i)*((F.F4(j-1,i) - F.F4(j,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G4(j-1,i) - G.G4(j,i))/Grid.delta_y_t));
        F_p.F1_p(j,i+1) = F.F1(j,i) + (dF_x.dF1_x(j)*delta_x);
        F_p.F2_p(j,i+1) = F.F2(j,i) + (dF_x.dF2_x(j)*delta_x);
        F_p.F3_p(j,i+1) = F.F3(j,i) + (dF_x.dF3_x(j)*delta_x);
        F_p.F4_p(j,i+1) = F.F4(j,i) + (dF_x.dF4_x(j)*delta_x);
    else % Forward difference with artificial viscosity
        dF_x.dF1_x(j) = (Grid.m1(j,i)*((F.F1(j,i) - F.F1(j+1,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G1(j,i) - G.G1(j+1,i))/Grid.delta_y_t));
        dF_x.dF2_x(j) = (Grid.m1(j,i)*((F.F2(j,i) - F.F2(j+1,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G2(j,i) - G.G2(j+1,i))/Grid.delta_y_t));
        dF_x.dF3_x(j) = (Grid.m1(j,i)*((F.F3(j,i) - F.F3(j+1,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G3(j,i) - G.G3(j+1,i))/Grid.delta_y_t));
        dF_x.dF4_x(j) = (Grid.m1(j,i)*((F.F4(j,i) - F.F4(j+1,i))/Grid.delta_y_t)) + (Grid.m2(j,i)*((G.G4(j,i) - G.G4(j+1,i))/Grid.delta_y_t));
        SF1 = ((0.6*(abs((P(j+1,i)) - (2*P(j,i)) + P(j-1,i))))/(P(j+1,i) + 2*P(j,i) + P(j-1,i)))*(F.F1(j+1,i) - (2*F.F1(j,i)) + F.F1(j-1,i)); % Value of artificial viscosity at j,i
        SF2 = ((0.6*(abs((P(j+1,i)) - (2*P(j,i)) + P(j-1,i))))/(P(j+1,i) + 2*P(j,i) + P(j-1,i)))*(F.F2(j+1,i) - (2*F.F2(j,i)) + F.F2(j-1,i));
        SF3 = ((0.6*(abs((P(j+1,i)) - (2*P(j,i)) + P(j-1,i))))/(P(j+1,i) + 2*P(j,i) + P(j-1,i)))*(F.F3(j+1,i) - (2*F.F3(j,i)) + F.F3(j-1,i));
        SF4 = ((0.6*(abs((P(j+1,i)) - (2*P(j,i)) + P(j-1,i))))/(P(j+1,i) + 2*P(j,i) + P(j-1,i)))*(F.F4(j+1,i) - (2*F.F4(j,i)) + F.F4(j-1,i));
        F_p.F1_p(j,i+1) = F.F1(j,i) + (dF_x.dF1_x(j)*delta_x) + SF1;
        F_p.F2_p(j,i+1) = F.F2(j,i) + (dF_x.dF2_x(j)*delta_x) + SF2;
        F_p.F3_p(j,i+1) = F.F3(j,i) + (dF_x.dF3_x(j)*delta_x) + SF3;
        F_p.F4_p(j,i+1) = F.F4(j,i) + (dF_x.dF4_x(j)*delta_x) + SF4;
    end
end
Updated_derivatives = dF_x;
Predicted_F_values = F_p;
end

So in this solution after expansion corner(x>E),the Mach no is maintained at more or less 2.2(2.18<M<2.2) at grid point 1(ie gridpoint located at wall where Abbot Boundary Condition is applied).But in my case( cfd code) one could see after expansion corner(x>E),Mach no at wall grid point keeps on decreasing,why it is happening,this is the problem in my code?And also only difference in my way is that i didnt predefine grid points along x direction and tried to follow it similar to a time marching procedure.What is the concept im lacking...