f[x0_,x1_,x2_,x3_]:=Sin[x1*x2+Cos[x0-x1]]+Sin[x2*x3+Cos[x1-x2]]


vec= D[f[x0,x1,x2,x3],{{x0,x1,x2,x3}}]
mat =D[vec, {{x0,x1,x2,x3}}]

{{-Cos[x0-x1] Cos[x1 x2+Cos[x0-x1]]-((Sin[x0-x1])^2) Sin[x1 x2+Cos[x0-x1]],Cos[x0-x1] Cos[x1 x2+Cos[x0-x1]]+Sin[x0-x1] (x2+Sin[x0-x1]) Sin[x1 x2+Cos[x0-x1]],x1 Sin[x0-x1] Sin[x1 x2+Cos[x0-x1]],0},{Cos[x0-x1] Cos[x1 x2+Cos[x0-x1]]+Sin[x0-x1] (x2+Sin[x0-x1]) Sin[x1 x2+Cos[x0-x1]],-Cos[x0-x1] Cos[x1 x2+Cos[x0-x1]]-Cos[x1-x2] Cos[x2 x3+Cos[x1-x2]]-(((x2+Sin[x0-x1]))^2) Sin[x1 x2+Cos[x0-x1]]-((Sin[x1-x2])^2) Sin[x2 x3+Cos[x1-x2]],Cos[x1 x2+Cos[x0-x1]]+Cos[x1-x2] Cos[x2 x3+Cos[x1-x2]]-x1 (x2+Sin[x0-x1]) Sin[x1 x2+Cos[x0-x1]]+Sin[x1-x2] (x3+Sin[x1-x2]) Sin[x2 x3+Cos[x1-x2]],x2 Sin[x1-x2] Sin[x2 x3+Cos[x1-x2]]},{x1 Sin[x0-x1] Sin[x1 x2+Cos[x0-x1]],Cos[x1 x2+Cos[x0-x1]]+Cos[x1-x2] Cos[x2 x3+Cos[x1-x2]]-x1 (x2+Sin[x0-x1]) Sin[x1 x2+Cos[x0-x1]]+Sin[x1-x2] (x3+Sin[x1-x2]) Sin[x2 x3+Cos[x1-x2]],-Cos[x1-x2] Cos[x2 x3+Cos[x1-x2]]-(x1^2) Sin[x1 x2+Cos[x0-x1]]-(((x3+Sin[x1-x2]))^2) Sin[x2 x3+Cos[x1-x2]],Cos[x2 x3+Cos[x1-x2]]-x2 (x3+Sin[x1-x2]) Sin[x2 x3+Cos[x1-x2]]},{0,x2 Sin[x1-x2] Sin[x2 x3+Cos[x1-xFortran2]],Cos[x2 x3+Cos[x1-x2]]-x2 (x3+Sin[x1-x2]) Sin[x2 x3+Cos[x1-x2]],-(x2^2) Sin[x2 x3+Cos[x1-x2]]}}

FortranForm[vec]

{-cos(x1 x2+Cos[x0-x1]] Sin[x0-x1],Cos[x1 x2+Cos[x0-x1]] (x2+Sin[x0-x1])-Cos[x2 x3+Cos[x1-x2]] Sin[x1-x2],x1 Cos[x1 x2+Cos[x0-x1]]+Cos[x2 x3+Cos[x1-x2]] (x3+Sin[x1-x2]),x2 Cos[x2 x3+Cos[x1-x2]]}