Algebraic System of Equations

Hello dear colleagues. I hope someone can help me.
I am solving a linear parabolic PDE on two variables, say, u1(x,t) and u2(x,t).
Then, for each node i want to calculate new quantities that depend on the values of u1 and u2. Those quantities are the solution of a nonlinear system of algebraic equations.
My code already solves perfectly de PDE itself but as i am not a good programmer, i don’t know how to write matrices and vectors of u1,u2 ( and its gradient) in C++ in order to perform a Newton method.
I want to do something like:

F = [f1(u1,u2) , f2(u1,u2)];
J = [ df1/du1, df1/du2; df2/du1, df2/du2];
dU = solve(F, -J);
u1 = u1 + du1;
u2 = u2 + du2;

Thank you very much.

Karl