How to plot elements from an array?

Question

Im trying to solve a differential equation for a range of parameters (1->4). I then want to automatically plot a graph for each of those parameters.

op = \[Rho]*Cp* D[T[t, r, z], t] -
   1/r*Subscript[k, th]*D[T[t, r, z], {r, 2}] -(*1/r*k*D[T[t,r,z],{r,
   1}]-*)
   Subscript[k, th]*(D[T[t, r, z], {z, 2}]) -
   sourcemulti[t, r, z, gap];

ic = {T[0, r, z] == 273};
(*the temperature of domain is Subscript[T, 0]=273.*)
\[CapitalOmega] = 
  Rectangle[{0, 0}, {20*beamradius, 
    10*1/\[Alpha]}]; (*define the domain of the solution*)
tend = 10*Subscript[\[Tau], 
   p]; (*simulation end time, 10*Subscript[\[Tau], p]*)

pde = {op == 0, ic};

array = Table[NDSolve[pde, {t, 0, tend}, {r, z} \[Element] \[CapitalOmega]], {gap, 1, 4}];

Table[Plot[array[[x]][t, 0, 1/\[Alpha]], {t, 0, tend}], {x, 1, 4}]

For each value of "gap" the NDsolveValue produces a function which is dependent on [t,r,z]. I then want to plot each of these functions. Ive tried combing it all under one "Table" function and a Do and For loop and none have worked so far. Do i need to make T a function of gap as well, or is their more to it? I know its only 4 graphs to plot but it will be expanded later. Any help appreciated im still learning mathematica

Edit. To make the code runnable you need a lof of constants and the source term. I tried to delete all of constants earlier to make it more streamlined but now ill dump the whole thing here:

\[Alpha] = 1.22*10^8;
    a = 0.0479;
    \[Rho] = 2700;
    Cp = 921;
    Subscript[\[Tau], p] = 1 10^-12; 
    n = 1.268 ;
    k = 10 ; 
    Subscript[k, th] = 236; (*thermal conductivity*)
    \[Lambda] = 1030 10^-9 ;(*wavelength in nm*)
    \[Alpha] = (4 \[Pi] k)/\[Lambda]; (*absorption coeff. m^-1*)
    R = ((n - 1)^2 + k^2)/((n + 1)^2  - k^2);(*reflectivity*)
    a = 1 - R; 
    beamradius = 
      10*10^-6 ;(*Beam has been focussed down to 10*10^-6 from 1.5*10^-9*)
    
    
    
    pulses = 2;
    pulse2start[gap_] := Subscript[\[Tau], p]*gap
    Subscript[I, 
      0] = (100 10^-6)/(Subscript[\[Tau], p] (10*10^-6/2)^2*\[Pi]);
    
    intensity[t_, r_] := 
     Piecewise[{{Subscript[I, 0], 0 <= r < beamradius}, {0, 
        beamradius <= r}}]
    Plot[intensity[t, 0], {t, 0, 2*Subscript[\[Tau], p]}]
    
    multipulse[t_, r_, gap_] := 
     intensity[t, r]* 
      Sum[Piecewise[{{0, t - (i - 1) pulse2start[gap] < 0}, {1, 
           t - (i - 1) pulse2start[gap] >= 0}}] - 
        Piecewise[{{0, 
           t - ((i - 1)*pulse2start[gap] + Subscript[\[Tau], p]) < 0}, {1,
            t - ((i - 1)*pulse2start[gap] + Subscript[\[Tau], p]) >= 
            0}}], {i, 0, pulses}]
    
    sourcemulti[t_, r_, z_, gap_] := 
     multipulse[t, r, gap]*a*\[Alpha]*Exp[-\[Alpha]*z]