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