Importing different plots onto one plot

I have a matlab code that prompts the user to enter specific details about an airfoil. then the user is able to plot different aspects of the airfoil. what i need to do is combine the plots from these two separate runs of the program. i have up loaded the entire code. I'm trying to explain it the best I can. basically, when you run the code for one type of airfoil, you get a certain result. when u run the code again for a different type of airfoil, new results are found, but the results from the initial airfoil are then gone. Is there a way to take the results from separate runs of the program and combine them into one plot? Thanks

clear all
clc

disp('Running MAE 424 vortex panel method code.');
disp(' ');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Part I: user input (airfoil geometry, # panels, min/max angles of attack)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

pmxx = input('Enter desired NACA 4-digit designation, e.g. 0012: ');
disp(' ');
Np = input('Enter number of vortex panels (even integer, ideally 40-100): ');
disp(' ');
alpha_min = input('Enter the minimum angle of attack to test (in degrees): ');
disp(' ');
alpha_max = input('Enter the maximum angle of attack to test (in degrees): ');
disp(' ');
kmax = input('Enter the max number of angles to test over this range (integer): ');
disp(' ');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Part II: creation of NACA 4-digit airfoil panely geometry
%
% Calculates the coordinates of the airfoil panel boundary points
% (xfoil,yfoil), and coordinates of the mean camber line points (xc,yc)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

disp('Creating (x,y) coordinates of airfoil boundary points, camber line...');
disp(' ');

pm    =  fix(pmxx/100);
xx    =  pmxx-pm*100;
p     = fix(pm/10);
m     = pm-p*10;
p     = p/100;
m     = m/10;
if m~=0
  fact1 = p/m^2; fact2 = 2*m; fact3 = p/(1-m)^2;
else
  fact1 = 0; fact2 = 19; fact3 = 0;
end
acoef = xx/20*[0.2969 -0.1260 -0.3516 0.2843 -0.1015];
%
%   Compute thickness distribution
%   vertices are clustered near the leading (x/c=0) and 
%   trailing (x/c=1) edges
%
beta = 0:2*pi/Np:2*pi; xc = 0.5*(1+cos(beta));
yt = acoef(1)*sqrt(xc)+xc.*(acoef(2)+xc.*(acoef(3)+...
        xc.*(acoef(4)+acoef(5)*xc)));
yt(Np/2+2:Np+1) = -yt(Np/2:-1:1);
%
%   Add the mean camber line
%
for i=1:Np+1
   if xc(i)<m
     yc(i) = fact1*xc(i)*(fact2-xc(i));
     tand  = fact1*(fact2-2*xc(i));
   else
     yc(i) = fact3*(1-fact2+xc(i)*(fact2-xc(i)));
     tand  = fact3*(fact2-2*xc(i));
   end
   cosd = 1/sqrt(1+tand*tand); sind = tand*cosd;
   xfoil(i) = xc(i)+yt(i)*sind; yfoil(i) = yc(i)-yt(i)*cosd;
end
% Plot airfoil panel geometry to check
figure(1)
plot(xfoil,yfoil,'r');
daspect([1 1 1]);
hold on
plot(xc,yc,'b');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Part III: calculations using vortex panel method
%
%    Calculates:  gamma - array of circulation densities;
%                 xcen  - ordinates of panel centers;
%                 cp    - pressure coefficient;
%                 clift - lift coefficient;
%                 cdrag - drag coefficient (aerodynamic);
%                 ista - starting panel i index;
%                 iend - ending panel i index;
%                 xcp - center of pressure;
%                 cmle - moment coeff. about leading edge;
%                 cmqc - moment coeff. about 1/4 chord.
%
%    All notations follow Kuethe & Chow (1998)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

disp('Running vortex panel method...');
disp(' ');

Np1 = Np+1;

alpha_min = alpha_min*pi/180;
alpha_max = alpha_max*pi/180;
alpha = linspace(alpha_min,alpha_max,kmax)'; % Angle of attack vector
alpha_deg = alpha*180/pi;

% Initialize column vectors of calculated quantities, fill with 0's
clift = zeros(kmax,1);
cdrag = zeros(kmax,1);
xcp = zeros(kmax,1);
cmle = zeros(kmax,1);
cmqc = zeros(kmax,1);
cl_exp = [.00395, .204, .308, .613, .917, 1.2065,1.4591];
a_exp = [.138,2.0496,2.9303,6.0206,8.8127,12.1,15];
for k = 1:kmax % Overall 'for' loop for varying angle of attack

%
%  Compute geometrical characteristics
%
for i=1:Np,
  ip1 = i+1;
        x(i) = 0.5*(xfoil(i)+xfoil(ip1));
        y(i) = 0.5*(yfoil(i)+yfoil(ip1));
        s(i) = sqrt((xfoil(i)-xfoil(ip1))^2+...
                    (yfoil(i)-yfoil(ip1))^2);
        theta(i) = atan2((yfoil(ip1)-yfoil(i)),...
                         (xfoil(ip1)-xfoil(i)));
        sine(i) = sin(theta(i)); cosine(i) = cos(theta(i));
        rhs(i) = sin(theta(i)-alpha(k));
end
%
%  Assemble the system matrix
%
for i=1:Np, for j=1:Np
   if (i==j)
     cn1(i,j)=-1;cn2(i,j)=1;ct1(i,j)=0.5*pi;ct2(i,j)=0.5*pi;
   else
     A =-(x(i)-xfoil(j))*cosine(j)-(y(i)-yfoil(j))*sine(j);
     B = (x(i)-xfoil(j))^2+(y(i)-yfoil(j))^2;
     C = sin(theta(i)-theta(j));
     D = cos(theta(i)-theta(j));
     E = (x(i)-xfoil(j))*sine(j)-(y(i)-yfoil(j))*cosine(j);
     F = log(1+s(j)*(s(j)+2*A)/B);
     G = atan2(E*s(j),B+A*s(j));
     P = (x(i)-xfoil(j))*sin(theta(i)-2*theta(j))+...
         (y(i)-yfoil(j))*cos(theta(i)-2*theta(j));
     Q = (x(i)-xfoil(j))*cos(theta(i)-2*theta(j))-...
         (y(i)-yfoil(j))*sin(theta(i)-2*theta(j));
     cn2(i,j) = D+0.5*Q*F/s(j)-(A*C+D*E)*G/s(j);
     cn1(i,j) = 0.5*D*F+C*G-cn2(i,j);
     ct2(i,j) = C+0.5*P*F/s(j)+(A*D-C*E)*G/s(j);
     ct1(i,j) = 0.5*C*F-D*G-ct2(i,j);
   end
end, end
%
%
%
for i=1:Np
   An(i,1)   = cn1(i,1); An(i,Np1) = cn2(i,Np);
   At(i,1)   = ct1(i,1); At(i,Np1) = ct2(i,Np);
   for j=2:Np
      An(i,j)=cn1(i,j)+cn2(i,j-1); At(i,j)=ct1(i,j)+ct2(i,j-1);
   end
end
An(Np1,1) = 1; An(Np1,Np1) = 1;
for j=2:Np, An(Np1,j) = 0; end
rhs(Np1) = 0; 
%
%  Solve for vortex panel strengths
%
gamma = An\rhs';
%
xcen = x;
i=1; while xcen(i+1)>0.99, i=i+1; end
ista = i;
i=Np; while xcen(i-1)>0.99, i=i-1; end
iend = i;
%
%  Compute Cp, Clift, and Cdrag
%
cx = 0; cy = 0; cmle(k) = 0; cmqc(k) = 0;
for i=ista:iend
   v(i) = cos(theta(i)-alpha(k));
   for j=1:Np1
      v(i) = v(i)+At(i,j)*gamma(j);
   end
   cp(i) = 1-v(i)^2;
   cx = cx+cp(i)*s(i)*sine(i);
   cy = cy-cp(i)*s(i)*cosine(i);
   cmle(k) = cmle(k)+cp(i)*s(i)*cosine(i)*xcen(i);
   cmqc(k) = cmqc(k)+cp(i)*s(i)*cosine(i)*(xcen(i)-0.25);
end
clift(k) = cy*cos(alpha(k))-cx*sin(alpha(k));
cdrag(k) = cy*sin(alpha(k))+cx*cos(alpha(k));
xcp(k) = -cmle(k)/clift(k);

end % for k = 1:

% Basic plot of lift coefficient vs. angle of attack
figure(2)
hold on
plot(alpha_deg,clift,'r')
xlabel('Alpha (deg)')

figure(2)
hold on
plot(alpha_deg,cdrag,'g')

figure(2)
hold on
plot(alpha_deg,cmqc,'b')

figure(2)
hold on
plot(alpha_deg,xcp,'m')

figure(2)
hold on
plot(a_exp,cl_exp,'k')

title('NACA 0012 - C_{drag},C_{lift},C_{m,1/4},X_{CP},C_{lift,exp} vs Alpha')
ylabel('C_{lift},C_{drag},C_{m,1/4},X_{cp},C_{lift,exp}')
legend('C_{lift}','C_{drag}','C_{m/1/4}','X_{cp}','C_{lift,exp}')

[per isakson]

There are many contributions in the File exchange, which address this problem, e.g. see Panel by Ben Mitch

Answer 1

[jridi Mariem]

commende subplot