I am slightly confused what help you would like from the Julia forum. I think the best thing would be to figure out why your Julia program produces discontinuous plots. The Julia solver routines, including the calls into LAPACK, seem to working correctly.
Perhaps to transfer numbers more smoothly between Julia and Fortran you should use ccall to call directly into your Fortran routine.
To help move this along and make it more transparent about what is happening, I made some modifications to your 21x21 code so that I wouldn’t have to do a bunch of commenting and uncommenting by just using deepcopy everywhere to be safe.
https://gist.github.com/mkitti/51cf831abf949d911de8220ca55da269
Also I modified the output as such:
@printf("Σ|x_0| = %g\n",sum(abs.(x0)))
@printf("Σ|x_1| = %g\n",sum(abs.(x1)))
@printf("Σ|x_2| = %g\n",sum(abs.(x2)))
@printf("Σ|x_3| = %g\n",sum(abs.(x3)))
@printf("Σ|x_4| = %g\n",sum(abs.(x4)))
@printf("Σ|x_5| = %g\n",sum(abs.(x5)))
for i=1:21
println(abs(x0[i]-x5[i]));
end
julia> include("21x21.jl")
Σ|x_0| = 9.34007e-06
Σ|x_1| = 6.80532e-19
Σ|x_2| = 2.99872e-18
Σ|x_3| = 5.97152e-19
Σ|x_4| = 1.97918e-19
Σ|x_5| = 7.5936e-19
1.8107496234338362e-19
6.214326071409102e-19
3.1071680867848654e-19
4.185954212237005e-19
9.354250834943223e-20
1.7618285302889447e-19
1.0369727249755868e-18
9.914084242815347e-19
7.589717713136635e-19
4.474386903525676e-19
1.9794114297145343e-18
1.0825916136606102e-18
4.642111499465803e-19
7.721670536909471e-19
2.9583648834882694e-19
1.205676349236012e-18
3.7738145890631644e-19
5.764605301733078e-19
5.611041877966055e-19
4.819371838949341e-19
9.340070946468615e-6