Sounds like in the linear case you got some vocab and denotations a bit confused...in G/FEM you don't discretize equations; what you have shown as "DISCRETIZATION:..." should just be calculus if you have done it right.
You are discretizing your solution using basis functions. You've indexed them as "j" within an element, so actually your elemental jacobian should be dRj/dci, unless by convention you want them flipped (which could be the case). So it might just be a matter of clarifying how you want everything to be indexed and keeping the indices straight?
In the linear case, you were not letting R be a matrix, that wouldn't make sense. You had no need for R at all. It would just have been the coefficients of your linear model, which would be teh same as dRj/dci if you had indexed everything correctly.
For each finite element, you would have a list (vector) of all c_j (the solution vector), and an associated residual vector, with number of entries equal to the number of basis functions (and therefore number of elemental nodes) NNE, and Jacobian is a matrix of size NNExNNE. At each iteration you are just solving the "linearized" version of your actual problem at the current solution vector guess.
Does that help? It's been a while since I've gone through all that myself, so maybe I'm not being super clear...
