How can I assign elements of an array in Modelica's equation and algorithm sections?

I am trying to solve a boundary value problem for two coupled PDEs by discritizing. For this I need to solve the two differential equations sequentially and iteratively loop until their values converge for a single timestep.

I thought to assign the boundary values in the equation section as they will always be valid. Because I need to iteratively solve the PDEs until the values converge, I implemented them in an algorithm section so they could be incorporated in a while loop.

However, due to Modelica specification 11.1.2 the entire array is initialized in the algorithm section and this causes the model to become singular and can't be solved, as there are too many equations (in this case there are 4 too many equations). Is there any way specify the boundary conditions for this problem, as this current method does not seem possible?

Here is a snippet of my code. I omitted the decleration of RHS variables for the statements in the algorithm section, but those would also be listed as satements. n and C are the variables relating to the PDEs I am trying to solve. They are arrays of length N, where N is the number of discritized elements.

equation 
  // PDE 1 Bounds 
  n[1] = n[2];  // dn/dx = 0 at X = 0
  n[N] = n[N - 1];// dn/dx = 0 at X = N

  // PDE 2 Bounds
    // Boundary Conditions
  C[1] = C[2];  // dC/dx = at X = 0
  C[N] = RH*Psat/(R*T);  // Boundary Condition at X = N. 

algorithm 
  //This section should run in a while loop   
    // PDE 1 Interior Node Calculation
  for i in 2:N-1 loop
    der(n[i]) :=(((D1[i + 1]*n[i + 1] + D1[i - 1]*n[i - 1] - 2*D1[i]*n[i])/dx^2)
       + PN*S[i]*Dv[i]*(aws[i]*Psat/(R*T) - C[i]))/(1 - V_g0);
  end for;

    // PDE 2 Exterior Node Calculation
  for i in 2:N-1 loop
    der(C[i]) :=(((Deff[i + 1]*C[i + 1] + Deff[i - 1]*C[i - 1] - 2*Deff[i]*C[i])
      /dx^2) + PN*S[i]*Dv[i]*(aws[i]*Psat/(R*T) - C[i]))/(V_g[i]);
  end for;

  //Update convergence criteria   

One possibility is to create helper arrays:

  Real smallN[N-2]=n[2:N-1];
  Real smallC[N-2]=C[2:N-1];
equation 
  // PDE 1 Bounds 
  n[1] = n[2];  // dn/dx = 0 at X = 0
  n[N] = n[N - 1];// dn/dx = 0 at X = N

  // PDE 2 Bounds
    // Boundary Conditions
  C[1] = C[2];  // dC/dx = at X = 0
  C[N] = RH*Psat/(R*T);  // Boundary Condition at X = N. 

algorithm 
  //This section should run in a while loop   
    // PDE 1 Interior Node Calculation
  for i in 2:N-1 loop
    der(smallN[i]) :=(((D1[i + 1]*n[i + 1] + D1[i - 1]*n[i - 1] - 2*D1[i]*n[i])/dx^2)
       + PN*S[i]*Dv[i]*(aws[i]*Psat/(R*T) - C[i]))/(1 - V_g0);
  end for;

    // PDE 2 Exterior Node Calculation
  for i in 2:N-1 loop
    der(smallC[i]) :=(((Deff[i + 1]*C[i + 1] + Deff[i - 1]*C[i - 1] - 2*Deff[i]*C[i])
      /dx^2) + PN*S[i]*Dv[i]*(aws[i]*Psat/(R*T) - C[i]))/(V_g[i]);
  end for;

  //Update convergence criteria   

But without the entire model I haven't checked it in detail.