So I have just started my journey learning Julia and really keen on Starting to use Gridap. I am still brushing up my knowledge on various topics.

Please refer to the below code, this code uses the differentialequations package to solve the 1D heat problem, very inelegantly. How could the same thing be reproduced using Gridap?

using DifferentialEquations
using Plots
function heat_equation(du,u,p,t)
# Parameters
alpha,dx = p
# Left boundary
du[1] = alpha*(-u[1]+u[2])/dx^2
# interior points
for i in 2:(length(u)-1)
du[i] = alpha *(u[i-1]-2*u[i]+u[i+1])/dx^2
end
# right boundary
du[end] = alpha*(u[end-1]-u[end])/dx^2
end
# u: tempreture T is a function of position x and time T
L = 1 #length of wire
dx = 0.01
x = 0.0:dx:L
# user inputs
left_temp = 0.9
right_temp = 0.1
alpha = 0.75
t_begin = 0
t_end = 1
# end user inputs
p = [alpha,dx]
tspan = (t_begin,t_end)
u_begin = [i < L/2 ? left_temp : right_temp for i in x]
prob = ODEProblem(heat_equation,u_begin,tspan,p)
sol = solve(prob)
plot(sol[:,34])

Note that finite-element methods are mainly about the spatial discretization of the problem. For a given spatial discretization, in principle there are many different time-stepping schemes you could use — you could even use DifferentialEquations.jl to do the time integration (also known as a method of lines).

(In your example code, you are also using the method of lines: combining a finite-difference spatial discretization with time integration via the ODE solver.)

Doing this well requires some understanding of both the spatial-discretization problem and the time-discretization problem. Heat equations are fairly stiff for time integration, so people often use an implicit scheme, like the Crank–Nicolson method (= theta method with θ=0.5) in the gridap example. DifferentialEquations has lots of other implicit schemes to choose from, and you can supply the matrix from Gridap as their Jacobian — you really want a sparse matrix here, so you don’t want to rely on standard autodiff.

(Of course, if you really have the heat equation with no time-dependent sources or materials, you can also just solve the time dependence with a matrix exponential.)