You need to save like shown bel and then plot:
clc
%% DESIGN ASSIGNMENT 1 %%
%%-----------given----------%%
Alt = 35000; %cruise altitude in ft
M = 0.8; %cruise mach
E = 1200; %loiter endurance in sec
Rd = 303805.774; %divert range in ft
Wc = 195*4; %total weight of crew(4) in lb
%%--------assumptions-------%%
Cc = 0.5/3600; %cruise SFC in per sec
Cl = 0.4/3600; %loiter SFC in per sec
b = 69.5; %wing span in ft
s = 520.4; %wing area in sq ft
AR = (b^2)/s; %calculated aspect ratio
S = 6.0; %ratio of s.wet to s
K = 13; %for high AR aircrafts
LDmax = K*AR/S; %maximum lift to drag ratio
LD = 0.866*LDmax; %lift to drag ratio
A = 1.02; %empty weight fraction constant
C = -0.06; %empty weight fraction constant
a = 969.16; %speed of sound in ft/sec at 35000 ft
V = M*a; %cruise velocity
%%-----------input----------%%
R = 5000000;
r = zeros(10,11) ;
w0 = r ;
wp = r ;
for i = 1:10
R = R + 500000 ;
x = 45;
r(i,j) = R ;
for j = 1:11
x = x + 5 ;
Wp = x*200; %total payload
wp(i,j) = Wp ;
%%---weight funtion ratio---%%
w1 = 0.97; %taxi-takeoff weight fraction
w2 = 0.985; %climb weight fraction
w3 = exp(-(R*Cc)/(V*LD)); %cruise weight fraction
w4 = 1; %decend weight fraction
w5 = exp(-(E*Cl)/LDmax); %loiter weight fraction
w6 = 0.995; %attempt to land weight fraction
w7 = 0.977; %climb weight fraction
w8 = exp(-(Rd*Cc)/(V*LD)); %divert cruise weight fraction
w9 = w4; %decend weight fraction
w10 = exp(-(E*Cl)/LDmax); %loiter weight fraction
w11 = 0.995; %landing weight fraction
wn = w1*w2*w3*w4*w5*w6*w7*w8*w9*w10*w11; %total weight fraction (w11/W0)
Wf = 1.06*(1-wn); %fuel weight fraction
%%-----take-off weight------%%
syms W0
eqn = W0 == (Wc + Wp)/(1-Wf-(A*((W0)^(C))));
gupta_suchita_DTOW = vpasolve(eqn,W0,1000000) ;
w0(i,j) = gupta_suchita_DTOW ;
end
end
mesh(w0,r)
mesh(w0,wp)
surf(w0,wp,r)
