Here is the full code with output and errors.

restart; with(plots): with(plottools):
w:=-Ps/2*(Hk*cos(theta)^2+2*H*cos(theta-psi));
             1    /             2                        \
      w := - - Ps \Hk cos(theta)  + 2 H cos(-theta + psi)/
             2                                            

dw:=diff(w,theta);
        1                                                         
dw := - - Ps (-2 Hk cos(theta) sin(theta) + 2 H sin(-theta + psi))
        2                                                         

equ:=dw=0;
                    1                                
           equ := - - Ps (-2 Hk cos(theta) sin(theta)
                    2                                

              + 2 H sin(-theta + psi)) = 0

Hg:=(sin(psi)^(2/3)+cos(psi)^(2/3))^(-3/2);
                                  1                  
           Hg := ------------------------------------
                                                (3/2)
                 /        (2/3)           (2/3)\     
                 \sin(psi)      + cos(psi)     /     

psi:=35.0*Pi/180;H:=0.1;Ps:=1.0;Hk:=1.0;
                      psi := 0.6108652383

                            H := 0.1

                           Ps := 1.0

                           Hk := 1.0

Hc:=evalf(Hk*Hg); theta:='theta';
                       Hc := 0.5103773339

                         theta := theta

theta_s:=evalf(Pi-arctan(tan(psi)^1/3));
                     theta_s := 2.912295114

H:='H';np:=100;H_max:=1.0;
                             H := H

                           np := 100

                          H_max := 1.0

h_p:=array(1..np):
h_n:=array(1..np):
m_p:=array(1..np):
m_n:=array(1..np):
h_p[1]:=0.0;m_p[1]:=cos(psi);
                          h_p[1] := 0.

                     m_p[1] := 0.8191520442

for i from 2 to np do
   h_p[i]:=(i-1)*H_max/(np-1):
   H:=h_p[i]:
   theta_plus:=fsolve(equ,theta,0..Pi/2):
   m_p[i]:=evalf(cos(theta_plus-psi)):
end do:

h_n[1]:=0.0; m_n[1]:=evalf(cos(Pi-psi));
                          h_n[1] := 0.

                    m_n[1] := -0.8191520446

for i from 2 to np-1 do
   h_n[i]:=(i-1)*Hc/(np-1):
   H:=h_n[i]:
   theta_minus:=fsolve(equ,theta,theta_s..Pi);
   m_n[i]:=evalf(cos(theta_minus-psi)):
end do:

h_n[np]:=Hc; m_n[np]:=evalf(cos(theta_minus-psi));
                    h_n[100] := 0.5103773339

    m_n[100] := cos(fsolve(1.000000000 cos(theta) sin(theta)

       + 0.5052220075 sin(theta - 0.6108652383) = 0, theta, 

      2.912295114 .. Pi) - 0.6108652383)

pl_pp:=array(1..np):
pl_pn:=array(1..np):
pl_np:=array(1..np):
pl_nn:=array(1..np):
i:='i';
                             i := i

NULL;

NULL;
NULL;
NULL;
for i from 1 to np do
   pl_pp[i]:=point([h_p[i],m_p[i]],color=red):
   pl_np[i]:=point([h_n[i],m_n[i]],color=blue):
   pl_pn[i]:=point([-h_n[i],-m_n[i]],color=red):
   pl_nn[i]:=point([-h_p[i],-m_p[i]],color=blue):
end do:
Error, (in plottools:-point) incorrect arguments for creating points structure, try providing the dimension option

display([seq(pl_pp[i],i=1..np),seq(pl_np[i],i=1..np),seq(pl_pn[i],i=1..np),seq(pl_nn[i],i=1..np)]);

Error, (in plots:-display) expecting plot structures but received: [pl_pp[60], pl_pp[61], pl_pp[62], pl_pp[63], pl_pp[64], pl_pp[65], pl_pp[66], pl_pp[67], pl_pp[68], pl_pp[69], pl_pp[70], pl_pp[71], pl_pp[72], pl_pp[73], pl_pp[74], pl_pp[75], pl_pp[76], pl_pp[77], pl_pp[78], pl_pp[79], pl_pp[80], pl_pp[81], pl_pp[82], pl_pp[83], pl_pp[84], pl_pp[85], pl_pp[86], pl_pp[87], pl_pp[88], pl_pp[89], pl_pp[90], pl_pp[91], pl_pp[92], pl_pp[93], pl_pp[94], pl_pp[95], pl_pp[96], pl_pp[97], pl_pp[98], pl_pp[99], pl_pp[100], pl_np[59], pl_np[60], pl_np[61], pl_np[62], pl_np[63], pl_np[64], pl_np[65], pl_np[66], pl_np[67], pl_np[68], pl_np[69], pl_np[70], pl_np[71], pl_np[72], pl_np[73], pl_np[74], pl_np[75], pl_np[76], pl_np[77], pl_np[78], pl_np[79], pl_np[80], pl_np[81], pl_np[82], pl_np[83], pl_np[84], pl_np[85], pl_np[86], pl_np[87], pl_np[88], pl_np[89], pl_np[90], pl_np[91], pl_np[92], pl_np[93], pl_np[94], pl_np[95], pl_np[96], pl_np[97], pl_np[98], pl_np[99], pl_np[100], pl_pn[59], pl_pn[60], pl_pn[61], pl_pn[62], pl_pn[63], pl_pn[64], pl_pn[65], pl_pn[66], pl_pn[67], pl_pn[68], pl_pn[69], pl_pn[70], pl_pn[71], pl_pn[72], pl_pn[73], pl_pn[74], pl_pn[75], pl_pn[76], pl_pn[77], pl_pn[78], pl_pn[79], pl_pn[80], pl_pn[81], pl_pn[82], pl_pn[83], pl_pn[84], pl_pn[85], pl_pn[86], pl_pn[87], pl_pn[88], pl_pn[89], pl_pn[90], pl_pn[91], pl_pn[92], pl_pn[93], pl_pn[94], pl_pn[95], pl_pn[96], pl_pn[97], pl_pn[98], pl_pn[99], pl_pn[100], pl_nn[59], pl_nn[60], pl_nn[61], pl_nn[62], pl_nn[63], pl_nn[64], pl_nn[65], pl_nn[66], pl_nn[67], pl_nn[68], pl_nn[69], pl_nn[70], pl_nn[71], pl_nn[72], pl_nn[73], pl_nn[74], pl_nn[75], pl_nn[76], pl_nn[77], pl_nn[78], pl_nn[79], pl_nn[80], pl_nn[81], pl_nn[82], pl_nn[83], pl_nn[84], pl_nn[85], pl_nn[86], pl_nn[87], pl_nn[88], pl_nn[89], pl_nn[90], pl_nn[91], pl_nn[92], pl_nn[93], pl_nn[94], pl_nn[95], pl_nn[96], pl_nn[97], pl_nn[98], pl_nn[99], pl_nn[100]]

NULL;

NULL;
 

 The following code:

for i to np do
    pl_pp[i] := point([h_p[i], m_p[i]], color = red);
    pl_np[i] := point([h_n[i], m_n[i]], color = blue);
    pl_pn[i] := point([-h_n[i], -m_n[i]], color = red);
    pl_nn[i] := point([-h_p[i], -m_p[i]], color = blue);
end do

Gives following error:

Error, (in plottools:-point) incorrect arguments for creating points structure, try providing the dimension option