Solving a large system of PDEs

Hello

I am currently working on trying to solve a large system of PDEs, more specifically, of Allen-Cahn type system, my approach is similar to the Cahn-Hilliard demo. So the system of PDE looks like this:

d n_i / dt = div( f_i(n_1, n_2, n_3, … n_k) ) for i = 1, 2, 3, … k
where f_i is a function that depends on all k of the n field

I have implemented it in FEniCS such that I could input k and it will solve the system of k equations. However, memory runs out pretty quick even when k is rather small (I could only do up to k=8 for my computer of 64GB RAM), and also the Newton solver iterates very slow too. I wander if there is any more efficient way to implement this solver than the approach in Cahn-Hilliard demo, this system will go as large as k=1000 so any advice on this will be really helpful, thanks!

(I tried breaking down the weak form L into k pieces and solve them individually, it wasn’t successful:( )