Just showing some RKs, no biggie.

rks.jl

const dz = 0.1
const z = 10
const space_steps = Int(z / dz)

const time_steps = 365 * 24 * 60
const ts = linspace(0, 12, time_steps)
const dt = (365 * 24 * 60 * 60) / time_steps

const alpha = 2e-6
const starting_temp = 30.0
const oscillations = 20.0

Ts = [starting_temp * ones(space_steps) for i in 1:time_steps]
for i in 1:length(Ts)
    Ts[i][1] = starting_temp + oscillations * sin.(1.0 * pi * ts[i] / 12)
end
const cache = similar(Ts[1])

function f!(out,T)
    adzinv = alpha * 1/dz^2
    @inbounds for i in 2:length(T)-1
        out[i] = adzinv * (T[i-1] - 2 * T[i] + T[i+1])
    end
end

function euler_solve!(Ts)
    for t in 1 : length(Ts)-1
        f!(cache,Ts[t])
        next_T = Ts[t+1]
        @inbounds for i in 2:length(next_T)-1
            next_T[i] = Ts[t][i] + dt * cache[i]
        end
        next_T[end] = Ts[t][end-1]
    end
end

using BenchmarkTools
@benchmark euler_solve!(Ts)

BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     141.518 ms (0.00% GC)
  median time:      158.612 ms (0.00% GC)
  mean time:        160.399 ms (0.00% GC)
  maximum time:     195.177 ms (0.00% GC)
  --------------
  samples:          32
  evals/sample:     1

const k1 = similar(cache)
const k2 = similar(cache)
const k3 = similar(cache)
const k4 = similar(cache)

function rk4_solve!(Ts)
    for t in 1 : length(Ts)-1
        T = Ts[t+1]
        Tprev = Ts[t]

        f!(k1,T)
        @inbounds for i in 2:length(T)-1
            T[i] = Tprev[i] + 0.5 * dt * k1[i]
        end

        f!(k2,T)
        @inbounds for i in 2:length(T)-1
            T[i] = Tprev[i] + 0.5 * dt * k2[i]
        end

        f!(k3,T)
        @inbounds for i in 2:length(T)-1
            T[i] = Tprev[i] + dt * k3[i]
        end

        f!(k4,T)

        @inbounds for i in 2:length(T)-1
            T[i] = Tprev[i] + dt * 1/6 * (k1[i] + 2k2[i] + 2k3[i] + k4[i])
        end
        T[end] = Tprev[end-1]
    end
end

using BenchmarkTools
@benchmark rk4_solve!(Ts)

BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     268.854 ms (0.00% GC)
  median time:      272.040 ms (0.00% GC)
  mean time:        272.894 ms (0.00% GC)
  maximum time:     281.205 ms (0.00% GC)
  --------------
  samples:          19
  evals/sample:     1


e_Ts = copy(Ts)
r_Ts = copy(Ts)


41.255129203639875
48.607491649368185

33.60947457082714
40.41290312408278
40.98357970304461

33.60947457082714
40.41290312408278
40.98357970304461
euler_solve!(e_Ts)
rk4_solve!(r_Ts)

println(Ts[100_001][1])
println(Ts[200_001][1])
println()
println(e_Ts[100_001][33])
println(e_Ts[300_001][48])
println(e_Ts[500_001][75])
println()
println(r_Ts[100_001][33])
println(r_Ts[300_001][48])
println(r_Ts[500_001][75])