Julia code to extract the coefficients of the heat equation.

poisson.jl

using ModelingToolkit

@parameters x y z
@variables u(..) # Need to choose u but user can pass the argument

Dx=Differential(x)
Dy=Differential(y)
Dz=Differential(z)

Dxx=Differential(x)^2
Dyy=Differential(y)^2
Dzz=Differential(z)^2

#DE
#eq=Dxx(u(x,y,z))+Dyy(u(x,y,z))+Dzz(u(x,y,z)) ~ 0
#eq=Dx(sin(x)*Dx(u(x,y,z)))+Dyy(u(x,y,z))+Dzz(u(x,y,z)) ~ 0
eq=Dx((x^2+exp(x)^2+y^2+sin(x)^2)*Dx(u(x,y,z)))+Dy(exp(x)*Dy(u(x,y,z)))+Dzz(3*u(x,y,z)) ~ 0

div_arg=Array{Any}(undef,length(eq.lhs.dict))
coeffs=Array{Any}(undef,length(eq.lhs.dict))
for m=1:length(eq.lhs.dict)
    # Extract the terms.
    T=collect(eq.lhs.dict)[m].first
    # Get the terms inside the divergence.
    # This is equivalent to throwing the derivative to a test function in the weak form
    div_arg[m]=(Int64(T.f==Differential(x))+Int64(T.f==Differential(y))+
                Int64(T.f==Differential(z)))*T.arguments[1]

    if(typeof(div_arg[m])==Term{Real,Nothing})
        #If the Differential has a constant function
        if(typeof(div_arg[m].arguments[1])==SymbolicUtils.Mul{Real,Int64,Dict{Any,Number},Nothing})
            coeffs[m]=Num(div_arg[m].arguments[1].coeff);
            continue
        end
        if(div_arg[m].arguments[1].f.name==:u)
            coeffs[m]=1.;
            continue
        end
    elseif(typeof(div_arg[m])==SymbolicUtils.Mul{Real,Int64,Dict{Any,Number},Nothing})
        # Else look for the function that multiplies the gradient.
        TERMS=collect(div_arg[m].dict)
        coeffs[m]=TERMS[2].first # Set coeffs to match type
        # Loop over terms
        for n=1:length(TERMS)
            # Check if there is a + sign in the term terms, then the terms are not separated
            if(typeof(TERMS[n].first)==SymbolicUtils.Add{Real,Int64,Dict{Any,Number},Nothing})
                coeffs[m]*=TERMS[n].first^(TERMS[n].second) # If yes, set the term as it is.
                break
            end
            # If product exists, terms are separated, so continue multiplying
            if(typeof(TERMS[n].first.f)!=Differential)
                coeffs[m]*=TERMS[n].first^(TERMS[n].second)
            end
        end
        coeffs[m]/=TERMS[2].first
        #print(coeffs[m],"\n")
    end
    count=1;
end

print("\nThe coefficients are: \n",coeffs)
print("\n\nThe corresponding divergence terms are: \n",div_arg,"\n\n")