int n = 0; //n is the dimensionality of our equations, nxn matrix
  float r[] = new float[N]; //residues vector b-Ax
  float x[] = new float[N]; //variables vector, will be returned into V1[] drawn as blue curve
  float p[] = new float[N]; //temp vector
  float q[] = new float[N]; //temp vector
  int xmap[] = new int[N]; //xmap represents the mapping from points on [0..N] range to [0..n] range. It is implicitly used for the matrix A.

void lorInit() {
 int index = 0;
  n=0; for(int i=0; i<N; i++) if (!C[i]) n++; //compute how many unknowns: sets n
  println("n="+n);
  for(int i=0; i<n; i++) p[i] = 0;  // step 2 pP=0
  for(int i=2; i<N-2; i++) // step 1: r=b-Ax
    if (!C[i]) {
      r[index] = -( 3*V1[i] -2*V1[i-1] -2*V1[i+1] +V1[i-2]/2 +V1[i+2]/2); // A/2-2B+3C-2D+E/2 
      xmap[i] = index;
      x[index] = V1[i];   // x is initialized to V and will evolve to the solution
      index++;            // next unknown      
      };
   }

boolean lorStep() {
    s = 0;
    for(int k=0; k<n; k++) s += r[k]*r[k];    
    if (s<0.000001) {println("s="+s+" too small"); stp='s'; return(false);   };        // early break (residues are zero)
    for(int k=0; k<n; k++) p[k] += r[k]/s;  // step 6 p=p+r/s
    for(int k=0; k<N; k++)    {
      if (!C[k])  {
        q[xmap[k]] = 3*p[xmap[k]];
        if (!C[k-2]) q[xmap[k]] += p[xmap[k-2]]/2;
        if (!C[k-1]) q[xmap[k]] += -2*p[xmap[k-1]];
        if (!C[k+1]) q[xmap[k]] += -2*p[xmap[k+1]];
        if (!C[k+2]) q[xmap[k]] += p[xmap[k+2]]/2;
      }
    }
    t = 0;  for(int k=0; k<n; k++) t += q[k]*p[k]; // step 8: this does tval = pTq
    if (s*t <= 6/1000) {println("s*t too small");  stp='t'; return(false);   };
    for(int k=0; k<n; k++) r[k] = r[k] - q[k]/t;  // step 10: r=r-q/t
    for(int k=0; k<n; k++) x[k] = x[k] + p[k]/t;   // step 11: x=x+p/t
    int index = 0; for(int i=0; i<N; i++)  if (!C[i]) V1[i] = x[index++];  // update V1
    index = 0; for(int i=0; i<N; i++)  if (!C[i]) R3[i] = r[index++];  // update V1
    stp=' '; return(true);
    }