Generating a plot from values

Question

I have the following function, which generates a value PO for a set of parameters w1,... w7. I would like to make a plot - change w4 alone and plot PO vs w4. How do I do this? In other words, how do I make a function out of this piece of code?

Remove["Global`*"]
li = 1000;
w1 = 0.04; w2 = 0.14; w3 = 0.1; w4 = 0.16; w5 = 0.08; w6 = 0.03; w7 = \
0.01; aji = 1; aki = 0.4; ami = 0.8;

w = Solve[{li*(w1 + w4*(1 - v) + w5*(1 - m) + w7*(1 - v)*(1 - m))*u ==
  aji, li*(w4*(1 - u) + w6*(1 - m) + w7*(1 - m)*(1 - u) + w2)*v ==
  aki, li*(w5*(1 - u) + w6*(1 - v) + w7*(1 - v)*(1 - u) + w3)*m ==
  ami}, {u, v, m}];
uu = w[[2, 1, 2]]; vv = w[[2, 2, 2]]; mm = w[[2, 3, 2]];
sol = Solve[{(w1 + w4*(1 - v) + w5*(1 - m) + w7*(1 - v)*(1 - m))*(u + 
        li*A) + 
     li*u*(-w4*B - w5*C - w7*(1 - v)*C - w7*(1 - m)*B) == 
    0, (w2 + w4*(1 - u) + w6*(1 - m) + w7*(1 - u)*(1 - m))*(v + 
        li*B) + 
     li*v*(-w4*A - w6*C - w7*(1 - u)*C - w7*(1 - m)*A) == 
    0, (w3 + w5*(1 - u) + w6*(1 - v) + w7*(1 - u)*(1 - v))*(m + 
        li*C) + 
     li*m*(-w5*A - w6*B - w7*(1 - u)*C - w7*(1 - m)*A) == 0}, {A, 
   B, C}] /. u -> uu /. v -> vv /. m -> mm;
A1 = sol[[1, 1, 2]]; B1 = sol[[1, 2, 2]]; C1 = sol[[1, 3, 2]]; 
DUDL = -(w1*uu + w2*vv + w3*mm + w4*(1 - (1 - vv) (1 - mm)) + 
 w5*(1 - (1 - mm)*(1 - uu)) + w6*(1 - (1 - vv)*(1 - uu)) + 
 w7*(1 - (1 - uu)*(1 - vv)*(1 - mm))) -  li*(w1*A1 + w2*B1 + w3*C1 + w4*(C1) (1 - vv) + w4*(B1) (1 - mm) + 
 w5*(C1) (1 - uu) + w5*(A1) (1 - mm) + w6*(A1) (1 - vv) + 
 w6*(B1) (1 - uu) + 
 w7*((1 - vv) (1 - mm) A1 + (1 - uu) (1 - vv) C1 + (1 - uu) (1 -mm) B1)) + aji/uu*A1 + aki/vv*B1 + ami/mm*C1
PO = -li*(w1*uu + w2*vv + w3*mm + w4*(1 - (1 - vv) (1 - mm)) + 
 w5*(1 - (1 - mm)*(1 - uu)) + w6*(1 - (1 - vv)*(1 - uu)) + 
 w7*(1 - (1 - uu)*(1 - vv)*(1 - mm))) + aji*Log[uu] +  aki*Log[vv] + ami*Log[mm]

Answer 1

Since solving the equations symbolically appears to take a very long time, instead I build a function that solves them numerically once a value for w4 is assigned:

ClearAll[po]
po[w4_?NumericQ] := 
 Module[{li, w1, w2, w3, w5, w6, w7, aji, aki, ami, w, uu, mm, sol, 
   A1, B1, C1, DUDL},
  li = 1000; w1 = 0.04; w2 = 0.14; w3 = 0.1;(*w4=0.16;*)
  w5 = 0.08; w6 = 0.03; w7 = 0.01;
  aji = 1; aki = 0.4; ami = 0.8;

  w = NSolve[{li*(w1 + w4*(1 - v) + w5*(1 - m) + w7*(1 - v)*(1 - m))*
       u == aji, 
     li*(w4*(1 - u) + w6*(1 - m) + w7*(1 - m)*(1 - u) + w2)*v == aki, 
     li*(w5*(1 - u) + w6*(1 - v) + w7*(1 - v)*(1 - u) + w3)*m == 
      ami}, {u, v, m}];

  uu = w[[2, 1, 2]]; vv = w[[2, 2, 2]]; mm = w[[2, 3, 2]];

  sol = NSolve[{(w1 + w4*(1 - v) + w5*(1 - m) + 
             w7*(1 - v)*(1 - m))*(u + li*A) + 
          li*u*(-w4*B - w5*C - w7*(1 - v)*C - w7*(1 - m)*B) == 
         0, (w2 + w4*(1 - u) + w6*(1 - m) + w7*(1 - u)*(1 - m))*(v + 
             li*B) + 
          li*v*(-w4*A - w6*C - w7*(1 - u)*C - w7*(1 - m)*A) == 
         0, (w3 + w5*(1 - u) + w6*(1 - v) + w7*(1 - u)*(1 - v))*(m + 
             li*C) + 
          li*m*(-w5*A - w6*B - w7*(1 - u)*C - w7*(1 - m)*A) == 0}, {A,
         B, C}] /. u -> uu /. v -> vv /. m -> mm;

  A1 = sol[[1, 1, 2]]; B1 = sol[[1, 2, 2]]; C1 = sol[[1, 3, 2]];

  DUDL = -(w1*uu + w2*vv + w3*mm + w4*(1 - (1 - vv) (1 - mm)) + 
       w5*(1 - (1 - mm)*(1 - uu)) + w6*(1 - (1 - vv)*(1 - uu)) + 
       w7*(1 - (1 - uu)*(1 - vv)*(1 - mm))) - 
    li*(w1*A1 + w2*B1 + w3*C1 + w4*(C1) (1 - vv) + w4*(B1) (1 - mm) + 
       w5*(C1) (1 - uu) + w5*(A1) (1 - mm) + w6*(A1) (1 - vv) + 
       w6*(B1) (1 - uu) + 
       w7*((1 - vv) (1 - mm) A1 + (1 - uu) (1 - vv) C1 + (1 - uu) (1 -
              mm) B1)) + aji/uu*A1 + aki/vv*B1 + ami/mm*C1;

  -li*(w1*uu + w2*vv + w3*mm + w4*(1 - (1 - vv) (1 - mm)) + 
      w5*(1 - (1 - mm)*(1 - uu)) + w6*(1 - (1 - vv)*(1 - uu)) + 
      w7*(1 - (1 - uu)*(1 - vv)*(1 - mm))) + aji*Log[uu] + 
   aki*Log[vv] + ami*Log[mm]
  ]

You can then plot po as a function of w4:

Plot[po[w4], {w4, -10, 0}]

Note that po seems to be real-valued only for $w4<0$, and complex-valued otherwise. If you are interested int he complex values, you will have to decide whether you want to plot their real and imaginary parts separately (look into ReIm for that), or their modulus, etc.