Graduate Diffusion equation and a system of equations with reciprocal unknowns?

[Hypatio]

TL;DR

Need to find an implicit numerical solution to a diffusion equation with concentration replaced with the derivative of a function of the concentration

So the normal diffusion equation looks like
\frac{\partial c}{\partial t} = k\frac{\partial}{\partial x}\left(\frac{\partial c}{\partial x}\right)
I know how to get explicit and implicit solutions to this equation using finite differences. However, I am trying to do the same for an equation of the following form:
\frac{\partial c}{\partial t} = k\frac{\partial}{\partial x}\left(\frac{\partial}{\partial x}\frac{\delta G}{\delta c}\right)
So, the term \delta G/ \delta c is a variational derivative of some function G. Spatial variation of this, instead of the concentration itself, acts as the force for diffusion. In my application, I need to find a solution in which G and \delta G/ \delta c are arbitrarily complicated, so I can probably settle for a solution in which it is a constant in time for each solve.

A 1D, forward in time, centered in space, finite-difference equation for an implicit solve might be:
c^1_i=c^0_i + \Delta t \frac{k}{\Delta x^2}\left(\phi^1_{i-1}-2\phi^1_{i}+\phi^1_{i+1}\right)
where superscript 0 indicates the old/reference time, and superscript 1 indicates the forward time being solved for, subscripts indicate position in space, and
\phi_i= \left(\partial G/ \partial c\right)_i
But maybe this is incorrect. I'm not sure how to write the equations so that the unknowns are the concentrations, c, and not phi, which is known. One approach I've tried is to write
\phi= \left(\delta G/ \delta c\right) = \frac{G-G_0}{c-c_0}= \frac{\Delta G}{c-c_0} =\frac{\Delta G}{c}
but this results in
c^1_i-\Delta t \frac{k}{\Delta x^2}\left(\Delta G_{i-1}\frac{1}{c^1_{i-1}}-2\Delta G_{i}\frac{1}{c^1_{i}}+\Delta G_{i+1}\frac{1}{c^1_{i+1}}\right) = c^0_i
which would result in a matrix of reciprocal and non-reciprocal unknowns, which I'm not sure how to solve for or rewrite to solve normally.

Any ideas about how to solve this type of equation implicitly?

Thanks for any help.

 

[Hypatio]

I've realized I may have made a mistake and I should have wrote the approximation of the derivative in a different way. For example
\phi_i = x_i+y_ic_i
for a first order approximation
\phi_i = x_i+y_ic_i+z_ic_i^2
for a second order approximation, or
\phi_i = x_i \ln (c)+y_i
for a logarithmic approximation.

because of the way c appears in the terms I can manage to write the coefficients that will go into the stiffness matrix. Does this make sense now?

On the other hand, I don't know how to deal with the squared term in the second order approximation, and even less how to deal with the logarithmic term in the third approximation. Any ideas?

 

[mathman]

This seems to be in the wrong forum. Should be in Differential Equations.